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1. INTRODUCTION 

Although from a historical point of view, composite materials have found practical use 
for centuries, during the past two decades there has been a tremendous increase in the use 
of composite materials in engineering applications, particularly in aerospace engineering. 
The main attraction of these materials is the ability to design and manufacture such 
materials to sustain a specific type of loading in the most efficient manner. If properly 
produced, composites can often achieve a combination of properties that are far superior 
to the properties of the individual constituents acting independently. 

It is well known that gas turbine engine structures, particularly those components di- 
rectly in the hot gas flow path, are subjected to extremely severe thermal and mechanical 
loading that can often lead to creep enhanced distortion, cracking and low cycle fatigue. 
As the demand for more efficient propulsion systems rise so does the thermal and mechan- 
ical loading. It is unlikely that the current generations of metal alloys would be suitable 
candidates for structural components in the future generations of efficient propulsion sys- 
tems. Ceramic components are often thought to be ideal as far as their thermal durability 
is concerned. Unfortunately, ceramics do not have adequate tensile strength to sustain a 
high level of mechanical loading. In recent years there has been significant effort in the 
attempt to incorporate fibrous inclusions within a ceramic matrix to develop a class of new 
materials, known as ceramic composites, for advanced engineering application. 

The mechanical behavior of ceramic composites under nonlinear, thermal, and dynamic 
loading is extremely complex and can only be understood if the observed behavior is 
interpreted in terms of micromechanical analyses. Such analyses must take care of the 
complex interaction of the individual fibers or bundles of fibers embedded in the three- 
dimensional ceramic matrix and must allow for increasing levels of sophistication in terms 
of the idealization of the fibers as well as the ceramic matrix. In addition complex interface 
behavior and controlled failure of the fiber must be considered. 
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It is evident that for proper micromechanical analysis of ceramic composites one needs 
to use a numerical method that is capable of idealizing the individual fibers or individual 
bundles of fibers embedded within a three-dimensional ceramic matrix. The analysis must 
be able to take account of high stress or temperature gradients from diffusion of stress or 
temperature from the fiber to the ceramic matrix and allow for the interaction between 
the fibers through the ceramic matrix. The analysis must be sophisticated enough to 
deal with failure of fibers described by a series of increasingly sophisticated constitutive 
models. Finally, the analysis must deal with micromechanical modeling of the composite 
under nonlinear thermal and dynamic loading. 

The boundary element method is uniquely suited for the task. BEM has proven its 
ability to accurately determine stress near a stress concentration. All functional quantities 
in a BEM system are on the boundary and interface surfaces, therefore, allowing nonlinear 
interaction on the interface between the matrix and the fiber to be readily described by fail- 
ure models. Furthermore, recent development has shown the generality and versatility of 
boundary element method in analyzing large two- and three-dimensional models subjected 
to static, dynamic, and thermal loads involving materials with nonlinear behavior. 

This report details progress made during the first four years towards the development 
of a boundary element code designed for the micromechanical studies of advanced ceramic 
composite. Additional effort has been made in generalizing the implementation to allow 
the program to be applicable to real problems in the aerospace industry. 

The ceramic composite formulations developed for this work have been implemented in 
the three-dimensional boundary element computer code ‘BEST3D’ which was developed 
for NASA by Pratt and Whitney and the State University of New York at Buffalo under 
contract NAS3-23697. BEST3D was adopted as the base for the ceramic composite pro- 
gram, so that many of the enhanced features of this general purpose boundary element 
code could be utilized. Some of these facilities include sophisticated numerical integration, 
the capability of local definition of boundary conditions, and the use of quadratic shape 
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functions for modeling geometry and field variables on the boundary. The multi-region 
implementation permits a body to be modeled in substructural parts; thus dramatically 
reducing the cost of the analysis. Furthermore, it allows a body consisting of regions of 
different ceramic matrices and inserts to be studied. 

In the next chapter the elastostatic and steady-state heat conduction boundary ele- 
ment formulation for ceramic composites is developed. The method for the semi-analytic 
integration of the kernel functions about the fiber is presented and the numerical imple- 
mentation of the formulation is discussed. This is followed by the development of the 
uncoupled thermoelastic BEM formulation in Chapter 3. Up to this point, only fibers 
assumed to be perfectly bonded to the composite matrix were considered. In Chapter 4, 
the previous formulations axe rederived in a form suitable for more general fiber to ma- 
trix interface connections. A spring and a nonlinear, coulomb friction constitutive models 
are developed for use as interface relations, and nonlinear solution algorithm is presented. 
Chapter 5 contains the transient heat conduction and transient uncoupled thermoelastic 
BEM formulations. In Chapter 6, nonlinear material behavior in the composite matrix is 
studied. Chapter 7 outlines the development of the general computer program implemen- 
tation. In Chapter 8, a number of numerical examples are presented to demonstrate the 
power of the present implementations. This report is then concluded with a summary and 
plan for future developments. 
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2. Steady-State Heat Conduction and Elastostatic BEM Formulation 

2.1 Introduction 

The conventional boundary integral equations for elastostatic and steady-state heat 
conduction analyses are used in deriving a boundary element formulation for the analysis 
of ceramic composite structures. The boundary integral equation written for a point £ in 
the interior of the composite matrix is modified by adding the boundary integral equations 
of each fiber written at the same point £ to this equation. This eventually eliminates the 
displacement (or temperature) variables on the fiber-matrix interface from the system, and 
therefore, reduces the total number of equations required for a solution of the system (see 
Figure 2-1). 

The BEM formulation for the steady-state heat conduction analysis of ceramic com- 
posites is identical to the elastostatic formulation, and therefore, the two formulations will 
be derived as one. 


2.2 Boundary Integral Equation Formulations 

The direct boundary integral equation for the displacement (or temperature) at a point 
£ inside an elastic composite matrix is 


= J [g£(x, £)<?(*) - F^(x,4)«f(x)]dS(x) 
+ E [<?£(*,0«f(*) - Fg(x,()u?(x) 


dS n {x) 


( 2 . 1 ) 


i,j =1,2,3 for elastostatics 
i,j = 1 for heat conduction 


where 

G™ , F™ are the fundamental solutions of the governing differential equations of the ceramic 
matrix of infinite extent 

Cij axe constants determined by the geometry at £ 
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Ui,u axe displacements and tractions (or temperature and flux) 

S, S” axe the surfaces of the outer boundary of the matrix and the n th hole (left for fiber), 
respectively 

N is the number of individual fibers 


Superscripts O and H identify the quantities on the outer surface of the matrix and the 
quantities on the surface of the hole, respectively. 

The conventional boundary integral equation for displacement (or temperature) can 
also be written for each of the N fibers. For the displacement (or temperature) at a point 
f inside the n th fiber we can write 


cfjiOMO = J sn [<?£(*, o*f(*) - f?(x,Ou?(x) 


dS"(x) 


( 2 . 2 ) 


i, j = 1,2,3 for elastostatics 
i,j = 1 for heat conduction 


GF,F£ are the fundamental solutions of the n th fiber 

*J V 

Cfj axe constants determined by the geometry at f in fiber n 
uf,tf are displacement and tractions (or temperature and flux) associated with the n th 
fiber 

S n the surface of the n th fiber 


Note, each fiber may have different material properties. 

We next examine the interface conditions between the composite matrix and the fiber. 
For a perfect bond the displacement (or temperature) of the matrix and the displacement 
(or temperature) of the fibers axe equal and the tractions (or fluxes) along the interface 
axe equal and opposite. 

«f(x) = wf(x) (2.3a) 

tf (*) = ~tf (x) (2.36) 

In elastostatics, when the elastic modulus of the fiber is much greater than the modulus 
of the composite matrix, the Poisson ratio of the fiber can be assumed equal to that of the 
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matrix with little error (no approximation is required for the heat conduction analysis). 
Therefore, upon consideration of the surface normals at the interface and examination of 
the Fij kernels, we can write the following relation for the n th fiber 


Fg(x,t) = -Fg(x,() 


(2.3 c) 


Substitution of equations (2.3) into equation (2.2) yields the following modified boundary 
integral equation for fiber n. 

CfjiOMt) = J s \~ G[j(x,Ot?(x) + ^(*,0«?(*)]dS"(x) (2.4) 

Finally adding the N fiber equations (2.4) to equation (2.1) and cancelling terms, yields 
the modified boundary integral equation for the composite matrix 


where 


C«(0m(0 = j s [g#(*,0<?(*) - F^(x, £)«?(*)] dS(x) 
+ E [ G«(x,tyt?(z)dS*(z) 

n=l Js " 


(2.5) 


5«(*,O = G0(*,O-G ! «(*.O 


Cy(^) are constants dependent on the geometry for a point f on the outer boundary and 
Cij(Z) = 6ij for a point £ in the interior of the body. 

2.3 Analytic Integration Around a Fiber 

The boundary element discretization of equations (2.4) and (2.5) in the conventional 
manner [Banerjee and Butterfield (1981)] requires a very fine discretization about the 
fiber/hole. Alternatively, a new formulation is introduced in this report for the efficient 
modeling and analysis of fibers/holes using what the authors refer to as ‘Fiber Elements’. 
The fibers are defined with Fiber Elements by describing the centerline of the (curvilin- 
ear, tubular) fiber with nodal points; defining the connectivity of the nodal points; and 
specifying the radius of the fiber at each of these nodal points (Figure 2.2). Internally 
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the program generates the surface of the fiber and the hole in which the resulting dis- 
placements (or temperatures) and tractions (or fluxes) axe described using a trigonometric 
circular shape function in the circumferential direction and a curvilinear shape function of 
any order in the longitudinal direction (the present work employs both linear and quadratic 
shape functions for this purpose). A long hole (which is allowed to vary in diameter) can 
be described by a number of the fiber elements connected end to end, and any fiber element 
not connected to another is assumed, by the program, to be closed at the end by a circular 
disc. 

Using the concept of the fiber element, the essential part of the formulation is the 
conversion of the two-dimensional surface integration of the fiber (and of the hole) to a 
one-dimensional integration. By performing a semi-analytical integration on the surface 
of the hole (or fiber) the numerical integration is significantly reduced. In equation (2.5) 
the integral under the summation is the integral which is to be modified. To facilitate an 
analytic integration in the circumferential direction, the three-dimensional kernel functions 
are first expressed in local coordinates with the center of the coordinate system coinciding 
with the center of the fiber/hole and the z axis aligned with the centerline of the fiber. 
The relative translation £< is added to the field coordinate & and the rotation is applied 
using the appropriate vector transformation. 

£»' = a ijij + 

where are the direction cosines between the axis of the local and global coordinate 
systems and the bar indicates a local variable. 

The integration point for a ring can now be expressed in cylindrical coordinates 
relative to the center of the fiber/hole as 

x\ — RcosO 

X2 = RsinB 

x 3 = 0 
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where R represents the radius of the fiber, i.e., R = (x 2 + xf) 1/2 . 

The normal vectors axe transformed by 

ni = n r cos9 

ri 2 = n r sin9 
n3 = n z 

where n r and n z represents the normals of the side of the hole in local coordinates and axe 
dependent on the change in the radius of the fiber/hole. On the side of a straight hole 
n r — 1 and n z = 0 , and on the flat surface closing the end of the hole/fiber n r = 0 and n z = 1. 

Next a circular shape function is employed to approximate the variation in the trac- 
tion (or flux) about the circumference of the fiber /hole. The circular shape function is 
multiplied and integrated with the three-dimensional Gij kernel, allowing the nodal values 
of traction (or flux) to be brought outside the integral. The shape function is expressed as 

ti = M^t] (summation over 7 is implied, 7 = 1, 2, 3) 

where 

M l (0) = ^ + jj cosB 
M 2 {6) = ^ -f sinO - ^ cosO 

O u O 

M 3 (6) = ^ p sinO - l-cosd 

w 3 3 3 

and <7 i s the nodal traction (or flux). See Figure 2.3. 

A modified circular shape function is used in the integration over the end of the hole to 

insure continuity of traction at the center of the end surface. The modified shape function 

is expressed as: 

M 7 = aM 7 + 6/3 7 = 1,2,3 

a = r/R b = (R — r)/R 

where 

R is the radius of the hole at the end, 
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AT 


is the location of the integration (Gauss) point as it sweeps from 
r = 0 to r = and 

is the circular shape function defined above. 


The traction must also be transformed between the local and the global systems by 


— Qiktk 


or 


tj — 

Since flux is a scalar quantity, no transformation is necessary in the heat conduction 
analysis. 

The last term in equation (2.5) can now be analytically integrated in the circumferential 
direction. For the m th hole the two integrals involved can be expressed as 

f G ij (x,Ot?(x)ds m (x)= f a jk f* G'£»\R,e,z,£)M^Rd9a ie t'ldC m 
Js m Jc rn Jo 

= [ G]j(R, z , £)tJdC m (z) 

Jc m 

where the indicated integration over C m is now a one- dimensional curvilinear integration 
along the hole and G? represent the analytically integrated fiber/hole kernels. Note, since 
the transformation vector is the independent of angle 9, it may be taken outside the d9 
integration. Similar analytic integration is also performed in equation (2.4). 

The final kernel functions of the fiber/hole obtained from the analytical integration 
are very lengthy and therefore are not presented in this report. They contain functions of 
elliptical integrals which in general are expressed numerically by common series approxi- 
mations. For a range of input values (coordinate locations), several higher order elliptical 
integral functions were found to produce incorrect numerical results. To overcome this 
problem, several new series were derived using a best fit polynomial approximation (as a 
function of the modulus of elliptic integrals) using values of the integrals calculated by a 
very accurate numerical integration in the circumferential direction. 
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The derivation of the fiber /hole kernels corresponding to the strain equation axe ac- 
complished from the displacement equation (2.5) by differentiation and application of the 
strain-displacement equations. The stress equation is then found using Hooke’s Law. Due 
to the complexity of the resulting fiber/hole kernels, the authors perform the required 
differentiation before the analytic circumferential integration. 

Finally we note, a fiber which has curvature along its length will differ in surface 
area about the circumference on the curved portion of the fiber. This is neglected in the 
formulation since the analytical integration is performed on an axisymmetric ring in which 
the surface area is constant about the circumference. This error, however, is small and 
disappears completely on a straight tubular fiber which is most commonly encountered. 
We should also mention that the fibers should not intersect the outer surface of the body 
or intersect other fibers. This minor restriction can be ignored if results at these locations 
axe not of interest. As long as the fibers do not coincide with nodal points of other elements 
the errors will be localized and will not affect the overall boundary element solution. 

2.4 Numerical Implementation 

The integral representations of Section 2.2 axe extremely accurate statements of the 
ceramic composite problem, however, approximations such as finite discretization and nu- 
merical integration axe necessary in order to obtain a solution to non-trivial problems. 
The goal of the numerical implementation of the present formulation is to obtain the most 
accurate and efficient implementation possible. 

2.4.1 Discretization 

After the analytical integration in the circumferential direction is complete, a fiber in a 
three-dimensional solid can be modeled as a two-dimensional curvilinear fine element with 
a prescribed radius at each longitudinal node. In the present work, linear and quadratic 
shape functions are utilized in modeling the geometry and field variables along the fiber 
element as well as the boundary elements on the outer surface of the 3-D body. In the 
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discretized form the displacement boundary integral equation for an elastic body containing 
fibers (equation (2.4)) can be expressed for a single fiber as 

p 


<%(()*( 0 = - E [ / G%(x,ON Q (r i)dCT 
P = 1 '■ Jcr 

+ E[ f Fi?(*,ON°(n)dC p 

p= 1 lJcf 


on f 




( 2 . 6 ) 


where 

P is the number of fiber elements, and 

N a (r]) represents a shape function over the curvilinear fiber element. Summation over 
a is implied. 

tf 7 and up are nodal values of traction (or flux) and displacement (or temperature) 
on the surface of the hole, respectively. 


In a similar manner, equation (2.5) can be discretized using one- and two-dimensional 
shape functions in the following manner. 

Q 


<?U(0«<(0 = E [/ s , G^(x,0L p (m,V2)dS^ 

-E 

+ E \f &lM,ON a {ri)dcAtT 

p= 1 


if 


( 2 . 7 ) 


where 

Q is the number of boundary elements on the outer surface of the composite matrix in 
the region, and 


lP(i h,rn) represents a two-dimensional shape function. Summation over /3 is implied. 


It is important to note that the displacement and traction (or temperature and flux) 
on a fiber /hole varies in the longitudinal as well as the circumferential direction, i.e., for 
displacement (or temperature). 

uj = M~ r N a u < p 
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The circular shape function M 7 has been analytically integrated into the kernel func- 
tions of equations (2.6) and (2.7). The ends of the fibers axe assumed to be a flat surface 
and a one-dimensional numerical integration is carried out in the radial direction. The 
coefficients obtained from the integration over the end are lumped with their respective 
coefficients from the integration of the side of the fiber. 

Equations for stress, strain, or flux at points in the interior of a body can be discretized 
and integrated in a similar manner. 

2.4.2 Numerical Integration 

The complexity of the integral in the discretized equation necessitates the use of nu- 
merical integration for their evaluation. The steps in the integration process for a given 
element is outlined below: 

1. Using appropriate Jacobian transformations, a curvilinear fiber element or boundary 
element is mapped on a unit line or on a flat unit cell, respectively. 

2. Depending on the proximity between the field point (£) and the element under con- 
sideration, there may be element subdivision and additional mapping for improved 
accuracy. 

3. Gaussian quadrature formulas are employed for the evaluation of the discretized inte- 
gral over each element (or sub-element). These formulas approximate the integral as a 
sum of weighted function values at designated points. The error in the approximation 
is dependent on the order of the (Gauss) points employed in the formula. To mini- 
mize error while at the same time maintaining computational efficiency, optimization 
schemes are used to choose the best number of points for a particular field point and 
element (Watson, 1979). 

4. When the field point coincides with a node of the element being integrated, the inte- 
gration becomes singular. In this case, the value of the coefficients of the Kj kernel 
corresponding to the singular node cannot be calculated accurately by numerical inte- 
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gration. Instead, after the integration of all elements is complete, this value is deter- 
mined so as to satisfy a rigid body displacement of the body (Banerjee and Butterfield, 
1981). 

2.4.3 Assembly of Equations 

After the derivation of the modified boundary integral equations and the analytical 
circumferential integration of the kernel functions, the next critical step in the formulation 
is the assembly of the fibers in the system equations. Here, efficiency is of utmost impor- 
tance. The approach to writing an efficient algorithm is to keep the number of system 
equations to a minimum by eliminating all unnecessary unknowns from the system. The 
strategy is to retain in the system only traction variables on the fiber-matrix interface. 
This is in contrast to a general multi-region problem where both displacement and trac- 
tions are retained on an interface. The elimination of the displacement on the interface is 
achieved through a backsubstitution of the fiber equations in the system equations which 
are made up exclusively from equations written for the composite matrix (on the outer 
surface and on the surface of the holes). The procedure is described below. 

Equation (2.7) is used to generate a system of equations for nodes on the outer surface 
of the composite matrix and for nodes on the surface of the holes containing the fibers. 
Written in matrix form we have 

On the Matrix Outer Surface: G M t° - F M u° + Gt H = 0 (2.8a) 

On the Matrix Hole Surface: G M t° - F M u° + Gt H = Iu H (2.86) 

where 

t° and u° are traction and displacement (or flux and temperature) vectors on the outer 
surface of the composite matrix 

t H and u H are traction and displacement (or flux and temperature) vectors on the hole 
I is the identity matrix 
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G m and F m matrices contain coefficients from the integration over the outer boundary. 
G matrix contains coefficients integrated about the fiber /hole 

Our goal is to eliminate u H from the system. To this end, equation (2.6) is written for 
every node on a fiber, collocating slightly outside the boundary of the fiber [at a distance 
of (1.25)*(fiber radius)] where C£(f) = 0. 

F F 2 u f = G Fa t F 

Superscript F2 identifies the equations written at points located slightly outside the bound- 
ary of the fibers. 

Noting u H = u F and t H = -t F we have 


F F 2 u h = — G F 2 t H 

Post multiplying equation (2.8b) by the F F2 matrix in equation (2.9) yields 

F F2 G M t° - F F 2 F m u° + F F2 Gt H = F F 2 u h 


( 2 . 9 ) 


( 2 . 10 ) 


Equation(2.9) can now be set equal to equation (2.10) and the final form of the system is 
derived. 

On Outer Surface: G M t° - F M u° + Gt H = 0 

On Hole: F F2 G M t° - F F2 F M u° + (F F2 G + G F2 )t H = 0 (2.11) 


At every point on the outer surface, either the traction or the displacement (or the 
temperature or the flux) is specified and on the surface of the hole only the tractions (or 
fluxes) are retained. Therefore, the number of equations in the system are equal to the 
final number of unknowns, and hence, the system can be solved. Thereafter, equation 
(2.8b) is used to determine the displacement on the fiber-matrix interface. 

It should be noted that since the displacement (or temperature) about a particular 
hole is present only in the fiber equation corresponding to that hole, backsubstitution can 
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be performed one fiber at a time in a more efficient manner than backsubstitution of all 
fibers at once. Further, nowhere in the assembly process is a matrix inversion necessary. 
This efficient assembly process was made possible due to the unique formulation of the 
modified boundary integral equations developed earlier in this chapter. 

When the composite matrix is divided into a multi-region model, the above fiber assem- 
bly is performed for each region independently. Thereafter, equilibrium and compatibility 
conditions are invoked at common interfaces of the substructured matrix composite. After 
collecting together the known and unknown boundary quantities, the final system can be 
expressed as 

A b x = B b y (2.12) 

where 

x is the vector of unknown variables at outer boundary unknown tractions (or fluxes) 
along the fiber-matrix interface 

y is the vector of known variables on the outer boundary of the composite matrix, 
and 

A b , B b are the coefficient matrices 

Standard numerical procedures are used to solve the unknowns in equation (2.12). 

■w 

Details are described in the computer development section (Chapter 6). 

Once the unknowns axe determined, the resulting displacements (or temperatures) 
on the fiber-matrix interfaces can be found using equation(2.8b) rendering a complete 
boundary solution for both the fibers and composite matrix. 

2.5 Interior Quantities 

Once all of the displacements and tractions (or temperature and fluxes) are known 
on the matrix outer surface and on the fiber-matrix interface, interior quantities of dis- 
placement, stress and strain (or temperature and flux) can be determined at any point in 
the composite matrix or fiber. For displacement (or temperature), either the conventional 
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boundary displacement (or temperature) integral equation (2.1) or (2.2) can be employed 
or alternatively the modified equations (2.3) or (2.5) can be used. 

Equations for strains can be derived from the forementioned displacement equations 
and the strain-displacement relations. Thereafter, equations for stress are obtained by 
substituting the resulting strain equations into Hooke’s law. 

The resulting equations, however, are not only invalid on the surface, but also difficult 
to evaluate numerically at points close to it. For points on the surface, the stresses can 
be calculated by constructing a local Cartesian coordinate system with the axes 1 and 
2 directed along the tangential directions and the axis 3 in the direction of the outward 
normal. The stresses ^ referred to these local axes (indicated by overbars) are then given 
by: 


v - Ev ( _ \ E 

011 = *3 + 2 ( £l1 + e 22 1 + T— : Eli 

1 — v \-v z \ ) 1 + 1/ 


012 — & 2 \ — 


E 


r £ i2 


2(1 + 1 /) 

v - Eu ( _ \ E 

022 = *3 + Z O £ll + £22 1 + r— £22 

i-i/ i-i/ z v / 1 + 1/ 

<*33 = <3 


<732 = 0 23 = ^ 2 
<731 = 013 = ti 


( 2 . 13 ) 


where E is the Young’s modulus, defines the components of the strains in the local axes 
system and are the traction on the boundary. This method of evaluating the stresses on 
the surface was originally devised by (Rizzo and Shippy, 1968). 
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II 








| a 
ji 
i 


\7 



= / 5n [G^(*,o<r(*)-^(« f o^(«)]ds"(«) 


Cy(0ui(0 = / s [g#(*>0*?(*) - Fg{x,t)u?{x) 
+ y~! f Gij(x,^)t l i { (x)dS ( x ) 

n^l 


dS(x) 


where 


G ij (x,0 = G? j (x,0-G( j (x,0 


i,j = 1,2,3 for elastostatics 
i,j = l for heat conduction 


Fig. 2.1 Boundary Integral Equation Formulation for Fiber Composite Materials 
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CIRCULAR SHAPE FUNCTION VALUE 


Fig. 2.2 The User Defines the Centerline of the Fiber and the Radius. The Surface of 
the Fiber (and the Surface of the Hole Containing the Fiber) is Automatically 
Generated by the Program 



Fig. 2.3 Value of the 3-Nodal Circular Shape Function about the Fiber/Hole 
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3. Steady-State Uncoupled Thermoelastic BEM Formulation 

3.1 Introduction 

The boundary element formulation for the steady-state uncoupled thermoelastic anal- 
ysis of perfectly bonded ceramic composite structures is similar to the elastostatic/heat 
conduction BEM formulation of Chapter 2. The primary difference is the inclusion of 
one-way coupling terms which, for the present formulation, requires the solution of the 
heat conduction equation prior to solving the thermoelastic equation set. This, however, 
is favorable since solving two uncoupled subsystems independently is more efficient than 
solving the system as a whole. 


3.2. Boundary Integral Equation Formulation 

Once again, the conventional boundary integral equation is the starting point for the 
uncoupled thermoelastic ceramic formulation. The displacement /temperature for a point 
£ inside the elastic composite matrix is 


Cij(0*W= I G%(z,t)t?(z)-F${z,t)u?( X ) 
Js 


dS(x) 


N 


+ 


u 




d5"(*) 


(3.1) 


i,j — 1, 2, 3, 4 


where 

G$f, Fjf are the uncoupled thermoelastic fundamental solutions of the governing differential 
equations of the ceramic matrix of infinite extent [Dargush, 1987] 

C,j are constants determined by the geometry at £ 

Ui,U are displacements and tractions for i = 1,2,3, and temperatures and fluxes for i = 4 
5, S n axe surfaces of the outer boundary of the matrix and the n th hold (left for fiber), 
respectively 

N is the number of individual fibers 
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Superscripts O and H identify the quantities on the outer surface of the matrix and the 
quantities on the surface of the hole, respectively. 

The conventional boundary integral equation for displacement /temper ature can also 
be written for each of the N fibers. For the displacement/temperature at a point £ inside 
the n th fiber we can write 

CfjiOMO = f gn - F£0r,O«f(*)]dS"(*) (3.2) 

hi = 1)2,3, 4 

Gfj, Fjj are the fundamental solutions of the n th fiber 

Cf 3 are constants determined by the geometry at £ in fiber n 
uf, if are displacement and tractions (i = 1 , 2, 3) and temperature and fluxes (i = 4) asso- 
ciated with the n th fiber 
S n the surface of the n th fiber 

Note, each fiber may have different material properties. 

We next examine the interface conditions between the composite matrix and the fiber. For 
a perfect bond the displacement /temperature of the matrix and the displacement /tempera- 
ture of the fibers along the interface are equal and the tractions/fluxes are equal and 
opposite. 

uf(x) = uf(x) (3.3a) 

t?(x) = -t[(z) (3.3 b) 

Substitution of equations (3.3) into equation (3.2) yields the following modified boundary 
integral equation for fiber n. 

CgiOMQ = J s \~ G ?A X ' 0*"(*) + *■"(*. 0«?(*) dS"(x) (3.4) 

Finally adding N fiber equations (3.4) to equation (3.1) and canceling terms, yields the 
modified boundary integral equation for the composite matrix 

• CijiOMO = J s [g£(*, 0<?(*) - ^(x,0«?(x)] dS(x) 
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dS^ix) 


(3.5) 


+ Xj J sn - Fii{x>Z)ui{x) 

where 

G «(*,0 = Go (*. 0 - G §(*,0 

Cy(^) are constants dependent on the geometry for a point ( on the outer boundary and 
Cij({) = Sij for a point in the interior of the body. 

3.3 Assembly of Equations 

The approach to writing an efficient algorithm is to keep the number of system equa- 
tions to a minimum by eliminating all unnecessary unknowns from the system. Once again 
the strategy used is to retain in the system only traction/flux variables on the fiber-matrix 
interface. This is in contrast to a general multi-region problem where both displace- 
ment/temperatures, and tractions /fluxes are retained on an interface. The elimination of 
the displacement/temperatures on the interface is achieved through a backsubstitution of 
the fiber equations into the system equations which axe made up exclusively from equations 
written for the composite matrix (on the outer surface and on the surface of the holes). 
The procedure is described below. 

Equations (3.4) and (3.5) can be discretized and integrated as described in chapter 2. 
The thermoelastic and heat conduction equations are integrated simultaneously. Therefore, 
using equation (3.4) a system of equations can be generated for the nodal points on the 
discretized composite matrix model. The equations for the nodes on the outer boundary 
of the matrix and for the nodes on the hole surface can be written separately in matrix 
notation as: 

On the Matrix Outer Surface: G M t° - F M u° + Gt H - Fu H = 0 (3.8a) 

On the Matrix Hole Surface: G M t° - F M u° + Gt H - Fu H = Iu H (3.86) 


(3.6) 

(3.7) 
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t° and u° are traction/flux and displacement /temperature vectors on the outer surface 
of the composite matrix 

t H and u H are traction /flux and displacement /temperature vectors on the hole 
I is the identity matrix 

G m and F M matrices contain coefficients from the integration over the outer boundary. 
G and F matrices contain coefficients integrated about the fiber/hole 

Our goal is to eliminate u H from the system. To this end, equation (3.4) is written for 
every node on a fiber, collocating slightly outside the boundary of the fiber [at a distance 
of (1.25)* (radius of fiber)] where C, F (f) = 0. 

G F2 t F - F F2 u f = 0 

Superscript F2 identifies the equations written at points located slightly outside the bound- 
ary of the fibers. 

Noting u H = u F , and t H = -t F we have 

F F2 u h = — G F2 t H (3.9) 

Post multiplying equation (3.8b) by the F F2 matrix in equation (3.9) yields 

F F2 G M t° - f F2 F m u° + F F2 Gt H - F F2 Fu h = F F2 u h (3.10) 

Equation (3.9) can now be set equal to equation (3.10) and the final form of the system is 
derived. 

On Outer Surface: G M t° - F M u° + Gt H - Fu H = O (3.11a) 

On Hole: F F2 G M t° - f F2 F m u° + (F F2 G + G F2 )t H - F F2 Fu H = O (3.116) 

Matrix F of equation (3.11) are coefficients derived from the integration of the h 
kernel defined in equation (3.7) as 

p. . — pff + f F 
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In the elastostatic/heat conduction formulation of Chapter 2 this kernel is zero according 
to equation (2.3c). This allowed the displacement /temperature vector u H to be removed 
from the formulation leaving an equal number of unknowns and equations for which a 
solution can be obtained. In the present formulation, if we once again assume the Poisson 
ratio of the fiber to be the same as the composite matrix, all but three of the sixteen 
components in the Fij kernel vanish. F 4 i,F 4 2 , and F 43 are non-zero when the coefficient of 
thermal expansion of the fiber is not the same as the matrix. Examining the Fy kernel 
more closely the following simplification is possible. 



■o 

0 

0 

Ft! 1 


a M -a F 

0 

0 

0 

F% 


a M 

0 

0 

0 

F% 



.0 

0 

0 

0 . 

UfJ 


(3.12) 


where 

are the displacement components of the hole 
axe the temperature component of the hole 
Fjtf, F™, Fj£ are components of the Fjf kernel on the hole surface of the matrix, and 

a F , a M are the coefficients of thermal expansion of the fiber and the matrix, respectively. 


From the above equation it is clear the displacement about the fiber /hole once again 
vanishes from the formulation, however, the temperature 6 H (0 H = u 4 ) is present as a 
coupling term, in the displacement equations (components 1,2 and 3 of the kernel) and 
vanishes in the temperature equation (component 4). Equation (3.11) may therefore be 


uncoupled and written as follows. 

Heat conduction system 

Outer Surface: G M t° - F M u° + Gt H = 0 (3.13a) 

Hole Surface: F F2 G M t° - F F2 F m u° + (F F2 G + G F2 )t H = 0 (3.136) 

Thermal elastic system 

Outer Surface: G M t° - F M u° + Gt H = B0 H (3.13c) 

Hole Surface: F F2 G M t° - F F2 F m u° + (F F2 G + G F2 )t H = B 2 0 H (3.13d) 
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where 


B0 h = Fu h 

and 

B 2 0 h = F F2 Fu h (3.146) 

In the heat conduction system the number of unknowns is equal to the number of equa- 
tions. The heat conduction system can be solved and the nodal temperatures on the hole 
determined using an uncoupled form of equation (3.8b). With the nodal temperatures on 
the hole known, the number of unknowns in the thermoelastic equation system is reduced 
to exactly the number of equations, and therefore, this system can be solved. With the 
use of equation (3.8b) all remaining nodal boundary variables can be determined on the 
fiber-matrix interface rendering a complete solution to the boundary value problem. Dis- 
placement, temperature, traction, flux, stress, or strain measures can now be determined 
at any point on the boundary or in the interior of the body as described in chapter 2. 
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4. Nonlinear Interface Connections Between the Fiber and the Matrix 

4.1 Introduction 

In order to accurately analyze a ceramic composite structure the interaction between 
the fiber and the composite matrix must be properly modeled. Interface phenomena such 
as perfect bonding, progressive debonding (with gap openings), frictional slipping, and 
plastic behavior along the fiber-matrix interface must be included. The failure mode of 
an interface is dependent on the state of stress at the interface. Therefore, the general 
mode of failure will be nonlinear and irreversible rendering a path dependent, quasistatic 
analysis. 

In this chapter the steady-state uncoupled thermoelastic boundary element formula- 
tions of chapter 3 is rederived in a form suitable for the inclusion of nonlinear interface 
connections. In addition to the perfectly bonded interface, two new types of interface 
connections are presented in a general form so as to be interchangeable upon input by the 
user. Finally the assembly of the numerically integrated BEM equations is presented and 
an incremental algorithm for their solution is described. The elastostatic and steady-state 
heat conduction formulations axe obtained from the uncoupled thermoelastic formulation. 


4.2. Boundary Integral Equation Formulation 

The conventional boundary integral equation for displacement/temperature is the 
starting point for the steady-state uncoupled thermoelastic ceramic formulation with non- 
linear fiber-matrix interface connections. The displacement/temperature for a point £ 
inside the elastic composite matrix is 


CW0«*(0 = J [<?#(*, 0*?(«) - ^(*,o«?(x)]ds(x) 

+ E / \g^(x, 0 ^(x)- F^ix^ixildS^ix) 

n= 1 Jsn I- -I 


( 4 . 1 ) 


i,j — 1 , 2 , 3,4 


where 
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G+j , F™ are the fundamental solutions of the governing differential equations of the ceramic 
matrix of infinite extent [Dargush, 1987] 

Cij are constants determined by the geometry at £ 

Ui,u are displacements and tractions 

5, 5” are the surfaces of the outer boundary of the matrix and the n th hole (left for fiber), 
respectively 

N is the number of individual fibers 


Superscripts O and H identify the quantities on the outer surface of the matrix and the 
quantities on the surface of the hole, respectively. 

The conventional boundary integral equation for displacement can also be written for 
each of the N fibers. For the displacement /temperature at a point £ inside the n th fiber we 
can write 


c£(0«i( 0= / 

Js n 


dSTix) 


(4.2) 


i,j = 1 , 2 , 3, 4 

GfjtFjj are the fundamental solutions based on material properties of the n th fiber 
Cfj are constants determined by the geometry at £ in fiber n 
uf , tf are displacement and tractions associated with the n th fiber 
S n the surface of the n th fiber 


Note, each fiber may have different material properties. 

We next examine the interface conditions between the composite matrix and the fiber. 
The difference between the displacement/temperature of the fiber and the displacement/ 
temperature of the fiber is defined by a variable di. The traction/flux on the hole is equal 
and opposite to the traction/flux on the fiber. Therefore 


uf (z) = uf (z) + di(x) 

(4.3a) 

<f (z) = -if (z) 

(4.36) 
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Substitution of equations (4.3b) into equation (4.2) yields the following modified boundary 
integral equation for fiber n. 

cS(0m(0 = J^[-GS(«, 0*^(0 -^(*, 0«f(*)]<*s"(x) (4.4) 

Finally adding N fiber equations (4.4) to equation (4.1) and invoking equation (4.3a), yields 
the modified boundary integral equation for the composite matrix 

N 

n=l 
N 

+ £ 


where 


= £ [<?£(*, #?(*) - /^(ar,0«f(x)]rfS(x) 

J sn ^ij(x,?)tf(x) - Flf(x,0di(x) - Fij(x,^)uf(x) dsr(x) (4.5) 


G ij (x,t) = G%(x,0-Gf j (x,0 
F ij (x,() = F t y(x,0 + FF(x,0 


For a point in the interior of the matrix 


Cij(0 = Sij and (C£(0)"=0 for all n 
For a point in the interior of the ifc th fiber 

C„(O = 0 ,nd (<$<{))” £ lit 

For a point on the outer boundary of the matrix 


Cij({) is dependent on the geometry at f, and 
((CfjiOT = 0 for all » 

For a point on the k th fiber-matrix 

c ij (o+((cr j (o) k = 6 i j 

((Gfj(0) n = 0 for n^k 


(4.6a) 


(4.66) 


(4.6c) 


(4.6d) 
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For a point £ on the k th matrix-fiber interface the left-hand side of equation (4.5) may be 
rewritten using equation (4.3a) and (4.6d) as 

cv(t)u?ao + (cs(0)\f(0 + £; (c$(o) Vo 

n=l 

njik 

= Cv(()uf(() + CijfOdifO + (c£«))\f(0 = + c«(0*(0 (4-7) 

4.3 Interface Constitutive Relationships 

4.3.1 Introduction 

The interface constitutive relationships of ceramic composite materials are, in general, 
nonlinear. This requires that their BEM solution be found using an incremental (qua- 
sistatic) algorithm. Therefore, the linear BEM equations of Section 4.2 are interpreted as 
incremental relations and the interface constitutive relationships are defined in terms of 
the incremental displacements/temperatures iij and the incremental tractions /fluxes u as 


= (4.8a) 

i,j = 1 to 3 for elastostatics 
i,j = 1 for heat conduction 
i, j = 1 to 4 for uncoupled thermoelasticity 


where k ^ is the nonlinear constitutive matrix dependent on the current state of stress, 


j • H 

dj — Uj uj , 


U = tf = ~i? 


(4.86) 

(4.8c) 


The elastostatic constitutive relationships that are derived in this section are expressed 
using a local coordinate system where the first component corresponds to the direction nor- 
mal to the interface and the second and third components correspond to arbitrary tangen- 
tial directions on the fiber-matrix interface. The local constitutive matrix is transformed 
to the global coordinate system when incorporated in the BEM equations as follows 


(k%f° hal 


= «*(*&) 


local 


Q mj 
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where a mj is the directional cosine transformation tensor. 


4.3.2 Linear Spring Interface Connection 

The general form of the linear spring interface connection may be expressed as (kfj 
replaces k^) 

U = k^dj ( 4 . 9 ) 

In the current elastostatic implementation the special case kf r = ib ni fcf 2 = jfc|g = k t , and kfj = 0 
for i ^ j, is assumed (see Figure 4 . 1 ). k n is the spring constant in the normal direction 
and k t is the spring constant in the tangential direction. The general form does not pose 
any more difficulty in implementation then the special case, but the additional generality 
is usually not required in practice. The limiting form as k n — *■ oo and k t —* oo approaches 
the case of a perfect bond, and k n = k t = 0 corresponds to a completely debonded interface. 
k n ~* o a and k t = 0 corresponds to a sliding interface with a perfect connection in the normal 
direction. Results using the present formulation with values k n = k t = (£ ,matrix + £ fiber ) * io 6 
compared very closely to results obtained for the same problem assuming a perfect bond. 
Similar comparisons were found to hold true for problems with jfc„ = k t = 0 verses E nber = 0 
(hole solution). 

In the heat conduction implementation k u = -k where k is the thermal conductivity, 
k = 0 does not permit heat flow across the interface which is analogous to an insulated hole 
and k — ► oo is analogous to a perfect bond between the fiber and the matrix. 

The kij for uncoupled thermoelastic analysis is derived by combining the elastostatic 
and heat conduction k^ to form a completely uncoupled k y matrix kn = k n , k 2 2 = £33 = 
k t , £44 = —k, and i fey = 0 for i ^ j. 


4.3.3 Spring-Friction Nonlinear Interface (Coulomb Friction) 

An interface model shown in Figure ( 4 . 2 ) exhibits a linear spring behavior normal to 
the interface when the normal tractions are in compression. Linear spring resistance is 
also observed in the tangential direction when the principal tangential traction is below 
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the slip limit defined by the coulomb friction criteria. When the slip limit is reached the 
traction that the interface can sustain in the tangential direction has reached a maximum 
and any additional tangential traction will result in an irrecoverable shift of the relative 
displacement between the fiber and the matrix, and a redistribution of stress (and therefore 
interface tractions) is required to bring the structure back in equilibrium. Furthermore, 
the model does not support tensile tractions normal to the interface. Instead the tractions 
at this point are set to zero through the constitutive relation kij = 0 (i,j = 1 to 3 ) and a 
gap between the fiber and the matrix will form. Once again a redistribution in stress is 
required to bring the structure back in equilibrium. 

The nonlinear interface constitutive relationship for a point on the interface exhibit- 
ing the spring-coulomb friction phenomena can be derived in a manner analogous to the 
incremental theory of plasticity [Selvadurai, 1988]. This requires a description and use of 

(1) A flow rule relating the irrecoverable nonlinear part of the displacement difference rate 
to the traction state at the interface. 

(2) A consistency relation requiring the new traction state to lie on the newly formed yield 
surface (defined by the hardening rule and yield function). 

(3) A yield function defining the limits of elastic behavior, and 

(4) A hardening rule defining the subsequent yield surfaces. 

The variable for the incremental displacement difference across the interface is assumed 
to be composed of an elastic part df and a plastic part d? 

di = d!-hdf (4.10) 

where the elastic component is related to the interface traction by the linear constitutive 
relation defined in Section 4.3.2. 

ii = kfjd* 
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or 


Next a flow rule is defined as 


i i = k? j (d j -d p j ) 




(4.11) 


3 . dtj 

where A is an unknown flow factor dependent on the current state of traction and Q is a 
nonlinear potential. Substituting the flow rule into equation (4.11) yields 


i — k e Ic e \ dQ 
** - k ij d J ~ 


(4.12) 


The consistency relation is defined as 


dF ■ „ 

= 0 

C/ti 


(4.13) 


Substitution of equation (4.12) into the consistency relationship and rearranging yields 
a relationship for A. 

1 AF 

(4.14) 


A — — — j k?d- 
Gdti '> 3 


where G = ^k e mn §^ 


Finally substitution of equation (4.14) into equation (4.12) yields 


ti = ki?dj 


(4.15) 


where 


and 


tep — ue _ hjp 
K ij ~ ij ij 


1 dF dQ 
^ G kiqkpj dt p dt q 


(4.16) 


(4.17) 


The nonlinear interface constitutive relationship is complete once F and Q are defined. 
Assuming the coulomb friction criteria, the yield function F is given by 


F = (t 2 2 + t 2 3 ) 1 ' 2 + nt n = 0 


(4.18) 
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where t n is the total traction normal to the interface, t 2 and t 3 are the total tangential 
tractions and n is the coefficient of friction between the fiber and the matrix. 

This relation holds only for t n < 0. If F < 0, the interface is assumed elastic and a 
spring connection is used. If F = 0 the nonlinear interface relation (4.15) is used. F < 0 
is meaningless in theory, however, since in practice we assume a finite size load step the 
function may become slightly positive in which case it is treated as F = 0. When t n > 0, 
kij = 0 (i,j = 1 to 3) is assumed. 

This relationship is assumed constant for the entire analysis and therefore to complete 
the analogy with classical plasticity, we can say we have zero hardening. The nonlinear 
potential Q is assumed to be the magnitude of the principal tangential traction. 

Q = (tl + il) 1/2 (4-19) 


If we assume the special elastostatic case described in Section 4.3.2 for kfj where kn = 
k n , k 22 = *33 = k t , and kfj = 0 for i ± j we can write for jfc?. 


j fcP. = 


0 

nk n t2 

l^krits 


0 

k t i 2 t 2 

k t t 2 t3 


0 

kti 2 iz 

ktisH 


where i 2 and t 3 axe the normalized components of total tangential tractions, i.e. 


(4.20) 


t7 ~ (ti+ky /2 and t3 ~ tt+k) 1 ' 2 

In the present heat conduction and uncoupled thermoelastic implementation, only a lin- 
ear (uncoupled) relation is allowed between the increment of flux and the temperature 
difference rate variable across the interface. This of course will be updated in the newly 
proposed work. 


4.4 Assembly of Equations for General Fiber-Matrix Interface Connections 

In the assembly of the perfectly bonded fibers the number of equations in the final 
system was reduced by eliminating all unnecessary unknowns from the problem. In the 
previous case all unknown displacement/temperature nodal variables on the fiber-matrix 
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interface were eliminated leaving only the unknown tractions/fluxes. In the present case the 
tractions/fluxes will be eliminated retaining only the fiber-matrix displacement difference 
variables along the interface. 

The steady-state uncoupled thermoelastic equations (4.4) and (4.5) can be discretized 
and integrated as described in Chapter 2. Therefore, using equations (4.4), (4.6c), and 
(4.7) a system of equations can be generated for the nodes on the discretized model of the 
composite matrix. The equations for the nodes on the outer boundary of the matrix and 
for the nodes on the hole surface can be written separately in matrix notation as: 

On the Matrix Outer Surface: G M t° - F M u° + Gt H - F H d - Fu f = 0 (4.21a) 

On the Matrix Hole Surface: G M t° - F M u° + Gt H - F H d - Fu f = Iu F (4.216) 

where 

t° and u° are traction/flux and displacement /temperature vectors on the outer surface 
of the composite matrix 

t H is the traction/flux vector on the hole 

u F is the displacement / temperature vector on the fiber 

d vector of displacement /temperature difference between the hole and fiber 

I is the identity matrix 

G M and F M matrices contain coefficients from the integration over the outer boundary 
G, F and F H matrices contain coefficients from the integration over the fiber /hole 

Vectors u° and d are arranged in the same order as the equations axe written so as to 
produce a strong main diagonal in matrices F M and F H . 

Equation (4.6c) is used in writing equation (4.21a). The Cy(£)ut(4) term is lumped with 
its respective (diagonal block) coefficient in matrix F M . The sum of these two terms axe 
accurately calculated indirectly using the rigid body technique [Banerjee and Butterfield, 
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1981]. Equation (4.7) is used in writing equation (4.21b). The Cij(£)di(£) term is lumped 
with its respective (diagonal block) coefficient in matrix F H (also calculated by the rigid 
body technique). Similarly, the term 6yuf(0 could be lumped with matrix F, but instead 
it is put on the right-hand side (forming an identity matrix I for use in subsequent matrix 
algebraic operations. Next, equation (4.5) is written for every shape function node on a 
fiber, collocating slightly outside the boundary of the fiber [at a distance of (1.25)*(radius 
of fiber)] where C£(£) = 0. 

G F2 t F - F F2 u f = 0 (4.22) 

Superscript F2 identifies the equations written at points located slightly outside the bound- 
ary of the fibers. Equation (4.22) can be rewritten as 

F F2 u f = — G F2 t H (4.23) 

Post multiplying equation (4.21b) by matrix F F2 of equation (4.23) yields 

F F2 G M t° - f F2 F m u° + F F2 Gt H - F F2 F H d - F P2 Fu f = F F2 u f (4.24) 

For efficiency, this operation can be carried out one fiber at a time. Equation (4.23) can 
now be set equal to equation (4.10). This yields, 

F F2 G M t° - f F2 F m u° + G 3 t H - F 3 d - F F2 Fu f = 0 (4.25) 

where 

G 3 = F F2 G -(- G fz , and 
F 3 = F F2 F h 

The nonlinear interface relation t H = -K ep d is substituted in equation (4.21a) and (4.25) 
and noting equations (3.12) and (3.14) the final form of the system can be written as 

On Outer Surface: G M t° - F M u° - (GK ep + F H )d = B 6 F (4.26a) 

On Hole Surface: F F2 G M t° - F F2 F m u° - (G 3 K ep + F 3 )d = B 2 0 F (4.266) 
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For each degree of freedom on every node on the outer surface, either the displace- 
ment/temperature or traction/flux is specified and on the fiber /matrix interface the dis- 
placement/temperature difference variable and the fiber temperature variable are present. 
Therefore the number of unknowns in the system is greater than the number of equations 
by exactly the number of nodes on the interface (exceeded due to the fiber temperature 
variables). However, as was the case in uncoupled thermoelasticity with a perfectly bonded 
fiber-matrix interface, the B and B 2 matrices influence only the thermoelastic solution, i.e., 
the coefficients corresponding to the heat conduction equations are zero. Therefore, in the 
heat conduction equations these terms do not appear. The heat conduction equation can 
be uncoupled from the system in a manner similar to equation (3.13) rendering an equal 
number of heat conduction equations and unknowns. The heat conduction system is solved 
and the fiber nodal temperatures are determined using equation (4.21b). With the fiber 
nodal temperatures known, the number of unknown in the thermoelastic equation system 
is reduced to exactly the number of equations, and therefore the system can be solved. 
With the use of equations (4.8) and (4.21b) all remaining nodal boundary variables can 
be determined on the fiber-matrix interface rendering a complete solution to the boundary 
value problem. Displacement, temperature, traction, flux, stress, or strain measures can 
now be determined at any point on the boundary or in the interior of the body as described 
in Chapter 2. 

The elastostatic and steady-state heat conduction formulations are derived from equa- 
tion (4.26). Noting the Fij kernel vanishes in elastostatic and heat conduction analysis, 
equation (4.21) is used setting matrix F = 0. The final system is similar to equation (4.26) 
with B = B 2 = 0. The number of unknowns are equal to the number of equations in a well 
posed problem and the system can be solved in one step. 

4.5 Numerical Algorithm for Nonlinear Interface Connections 

The nonlinear constitutive relationships are incorporated in the linear BEM equations 
and solved using the incremental (quasistatic) algorithm described below. 
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Nonlinear Solution Algorithm 

(a) The boundary (and body force loading if present) is divided into a number of small 
sub-increments. Ten sub-increments have been found to be sufficient. 

Initial Load Step (Iterative Start-up Loop): 

(b) Assemble the BEM system equation (4.26) as described in Section 4.4 using the in- 
terface relations calculated in step (d). In the initial assembly step, linear spring 
connections (Section 4.3.2) are assumed at each interface node. 

(c) Apply the sub-increment of load to the system equations (4.26). 

(d) Determine the traction state at each interface node and derive the appropriate interface 
constitutive relation for that node. 

- If the normal traction a the interface is positive, assume a gap opening with zero 

traction, i.e., = 0 . 

- If F < 0 assume a spring connection (Section 4.3.2) 

- If F < 0 assume a nonlinear connection (Section 4.3.3) and evaluate the constitutive 
relationship using the current traction state (from current loading). 

(e) Changing the interface connections from a linear spring to nonlinear relationship in the 
first sub-load step may cause a major redistribution in the tractions on the interface. 
Therefore, the constitutive relations calculated in step (c) will change. Repeat step 

(b) , (c) and (d) until a coverage solution is achieved for the first load step. 

Subsequent Load Steps: . 

(f) Steps (b), (c), and (d) are repeated once for each additional load step. In general, 
iteration within a subsequent sub-load step is not required since the constitutive re- 
lationships at interface nodes change in a gradual manner once the start-up loop is 
complete. The solution of any sub-increment is the accumulation of all previous sub- 
load steps. It is this current solution that is used in step (d) to evaluate the current 
interface constitutive relationships used in the next sub-incremental load step (step 

(c) )- 
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5. TRANSIENT HEAT CONDUCTION AND TRANSIENT UNCOUPLED 
THERMOELASTIC BEM FORMULATIONS 

5.1 Introduction 

The analysis of a body in a transient state is inherently more complicated than a 
steady-state analysis. The boundary element formulation for transient analysis contains 
convoluted integrals which must be integrated in time as well as in space. Therefore, a 
transient BEM analysis is expected to be more complicated and more expensive than its 
steady-state counterpart. To compound these difficulties, the cost saving techniques devel- 
oped for the steady-state ceramic composite analyses cannot be employed for the transient 
case. First, the efficient modified boundary integral formulations and the resultant as- 
sembly schemes of Chapters 2-4 do not lend themselves to transient fiber formulations 
for composites with general material properties. Secondly, due to the complexity of the 
transient kernel functions, the integration about the circumference of the fiber must now 
be carried out using numerical integration, thus adding to the cost of the analysis. 

Nevertheless, an efficient numerical integration scheme as well as an efficient assembly 
algorithm have been derived for the transient heat conduction and thermoelastic BEM 
analyses. In the assembly algorithm the equations of the fibers are backsubstituted into 
the equation system of the composite matrix, hence, reducing the size of the overall system. 
This is carried out prior to decomposition at a point when the equations for each fiber 
can be handled individually in the most efficient manner. This assembly, as well as the 
decomposition, is required only for the initial time step. Later time steps require only a 
new calculation of the system’s right-hand side. 

5.2 Transient Boundary Integral Equation Formulation 

The transient, uncoupled thermoelastic boundary integral equation (Dargush and Baner- 
jee, 1991) for the displacement (and temperature) at a point £ in a composite matrix is 
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( 5 . 1 ) 


Cpc(0 Mt> *) = J S {9 p a * u°{x, <)] ds(x) 

M . 

+ £ / [9&*tZ(X,t)-f£*u%(X,t)}dS"{X) 

m = 1 Jsm 


where 

a, f3 indices varying from 1 to 4 for uncoupled thermoelasticity and equaled to 4 for 
heat conduction 

Cpo, constants determined by the geometry at £ 
up,tp generalized displacement and traction 
Up = [til, u 2 , u 3 , 9] T 
tp = [< 1 , < 2 , * 3 ! q] T 

9, q temperature, heat flux 

9ap> fap generalized displacement and traction kernels 
S,S m are the surfaces of the outer boundary of the matrix and the mth hole (left for 
fiber) 

M is the number of individual fibers 
and, for example 

9pa *tp = f 9p* a (X,*',Z,T)tp(X,T)dT 
JO 

denotes a Riemann convolution integral. 

Superscripts O and H axe used to highlight the quantities associated with the outer 
surface of the matrix and the quantities associated with the surface of the hole, respectively. 

The boundary integral equation for mth fiber can be written as 


CpaiOMO = I [spa * - fp a * «£(*,*)] dS"(X) ( 5 . 2 ) 

JS m 

where 

9pa’ fp a are the fundamental solutions for a body with material properties corresponding 
to the mth fiber. 
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Cp a are constants determined by the geometry at £ in fiber m 
Up, ip are the generalized displacement and tractions associated with the mth fiber 

5 m the surface of the mth fiber 

In principle, at each instant of time progressing from time zero, these equations can be 
written at every point on the respective boundary. The collection of the resulting equations 
could then be solved simultaneously, producing exact values for all the unknown boundary 
quantities. In reality, of course, discretization is needed to limit this process to a finite 
number of equations and unknowns. Techniques useful for the discretization of (5.1) and 
(5.2) axe the subject of the following section. 

5.3 Numerical Implementation 

The boundary integral equations (5.1) and (5.2) developed in the last section, are exact 
statements. No approximations have been introduced other than those used to formulate 
the boundary value problem. However, in order to apply (5.1) and (5.2) for the solution 
of practical engineering problems, approximations are required in both time and space. In 
this section, an overview of a general-purpose, state-of-the-art numerical implementation 
is presented for transient analysis. 

5.3.1 Temporal Discretization 

Consider, first, the time integrals represented in (5.1) and (5.2) as convolutions. Clearly, 
without any loss of precision, the time interval from zero to t can be divided into N equal 
increments of duration at At. 

l 

By assuming that the primary field variables, tp and up , are constant within each At 
time increment, these quantities can be brought outside of the time integral. That is, 

N rn&t 

9pa*tp(X,t) = 'y't%(X) gpa(X - - r)dr (5.3a) 

N rn&t 

fp a *up(X,t) = T / u^(X) fp Q (X-U-r)dr (5.36) 

£rl J(n- l)At 
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where the superscript on the generalized tractions and displacements represents the time 
increment number. Notice, also, that within an increment, these primary field variables 
are now only functions of position. Since the integrands remaining in (5.3) are known in 
explicit form from the fundamental solutions, the required temporal integration can be 
performed analytically, and written as 

rnAt 

G$+ 1 - n (X-t)= / gfic(X-t,t-T)dT (5.4a) 

./(n-l)At 

F N+i-n {x _ 0= / f Pa (X-U-r)dr. (5.46) 

J(n— l)At 


These kernel functions, Gp a (X -£) and Fg a (X -£), are detailed in Appendix B. Combining 
(5.3) and (5.4) with (5.1) and (5.2) produces 

<W0«#( 0 = £ { / [< +1 ~ n (x - - F " a +1 ~ n (x - 0«3(*)] dS ( x ) 

n=l ^ Js 

M. \ 

+ £ m [C' n (^W-C' n ( X -W]^W (5.5a) 

m=l J S m ' 


and for the mth fiber 

<wo«*(o = £ { J sm K* +1 ' n (* - twm - F " a +l - n ( x - o«2(*)] ^ w} ( 5 - 5fc ) 

which are the boundary integral statements after the application of the temporal discretiza- 
tion. 


5.3.2 Spatial Discretization 

With the use of generalized primary variables and the incorporation of a piecewise 
constant time stepping algorithm, the boundary integral equation (5.5) begins to show a 
strong resemblance to the equations of the steady-state analyses, particularly for the initial 
time step (i.e., N = l). In this subsection, those similarities will be exploited to develop 
the spatial discretization for the uncoupled quasistatic problem. This approximate spatial 
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representation will, subsequently, permit numerical evaluation of the surface integrals ap- 
pearing in (5.5). The techniques described here, actually, originated in the finite element 
literature, but were later applied to boundary elements by Lachat and Watson (1976). 

The process begins by subdividing the entire surface of the body (including the fibers) 
into individual elements of relatively simple shapes. The geometry of each element is, then, 
completely defined by the coordinates of the nodal points and associated interpolation 
functions, as described in Chapter 2. 

In the present work, the geometry is exclusively defined by quadratic shape functions. 
On the other hand, the variation of the primary quantities can be described, within an ele- 
ment, by either quadratic or linear shape functions. (The introduction of linear variations 
proves computationally advantageous in some instances.) 

Once this spatial discretization has been accomplished and the body has been subdi- 
vided into Q boundary elements and P fiber elements, the boundary integral equation can 
be rewritten for the matrix as 

CpaiO^i 0 = E { E 
«= 1 9=1 


J SI 


-iN+l—n 


(xiO-OUOds^xic)) 


i n 


and 


- E [J st F" + 1 -"(X (0 - OMOrfsWO) 

+ E \J SP Gpa'^iXiO - t)NM)^( 6 )dS?(X( 0 ,e)\ 

-m. F^- n (x(0-0Nu,(0My(e)d^x(0,e)^ 


Cflaityrfit) = E { E f / G pa l ~ n (m - oauo^ws'WO, *)] t%^ 

n=l *■ p= 1 '-•'S? J 

- E [J SP F 0 * +l - n ( X ’ (C) - WUQM-rWdSPiXiC), 0)] ( 5 . 6 b) 


(5.6a) 


where 


Q is the number of boundary elements on the outer surface of the composite matrix 
in the region, 

P is the number of fiber elements in the region, 

L u represents a two-dimensional shape function, 

N 1 represents a one- dimensional shape function, 

M 7 is the circular shape function defined in Chapter 2, 

t, u are nodal values of generalized traction and displacement, respectively. 

In the above equation, the nodal quantities and u pw were brought outside the surface 
integrals. This positioning of the nodal primary variables outside the integrals is, of course, 
a key step since now the integrands contain only known functions. However, before dis- 
cussing the techniques used to numerically evaluate these integrals, a brief discussion of 
the singularities present in the kernels G pa and F po is in order. 

The fundamental solutions to the uncoupled quasistatic problem contain singularities 
when the load point and field point coincide, that is, is when r = 0. The same is true of G^ a 
and Fg a , since these kernels are derived directly from the fundamental solutions. Series 
expansions of terms present in the evolution functions can be used to deduce the level of 
singularities existing in the kernels. 

A number of observations concerning the results of these expansions should be men- 
tioned. First, as would be expected has a stronger level of singularity than does the 
corresponding G l a/3 , since an additional derivative is involved in obtaining F* p from G l ap . 
Second, the coupling terms do not have as high degree of singularity as do the correspond- 
ing non-coupling terms. Third, all of the kernel functions for the first time step could 
actually be rewritten as a sum of steady-state and transient components. That is, 

G 1 88 s-i , tr/~i 1 

a(3 ~ ° ap + ^a/3 

T?1 88 T? i tV Z7*l 

?aP — + * a(3 

Then, the singularity is completely contained in the steady-state portion. Furthermore, 
the singularity in Gjj and is precisely equal to that for elastostatics, while G l ee and Fg e 
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singularities axe identical to those for potential flow. This observation is critical in the 
numerical integration of the F a p kernel to be discussed in the next subsection. However, 
from a physical standpoint, this means that, at any time t, the closer one moves toward the 
load point, the closer the quasistatic response field corresponds with a steady-state field. 
Eventually, when the sampling and load points coincide, the quasistatic and steady-state 
responses are indistinguishable. As a final item, after careful examination of Appendix 
B, it is evident that the steady-state components in the kernels G^p and F£p, with n > 1, 
vanish. In that case, all that remains is a transient portion that contains no singularities. 
Thus, all singularities reside in the aa G a p and 88 F a p components of G l a p and F X p, respectively. 

5.3.3 Numerical Integration 

Having clarified the potential singularities present in the coupled kernels, it is now 
possible to consider the evaluation of the integrals in equation (5.6). 

To assist in this endeavor, the following three distinct categories can be identified. 

(1) The point £ does not lie on the element m. 

(2) The point £ lies on the element m, but only non-singular or weakly singular integrals 
are involved. 

(3) The point lies on the element m, and the integral is strongly singular. 

In practical problems involving many elements, it is evident that most of the integration 
occurring in equation (5.6) will be of the category (1) variety. In this case, the integrand 
is always non-singular, and standard Gaussian quadrature formulas can be employed. So- 
phisticated error control routines are needed, however, to minimize the computational 
effort for a certain level of accuracy. This non-singular integration is the most expensive 
part of a boundary element analysis, and, consequently, must be optimized to achieve an 
efficient solution. In the present implementation, error estimates, based upon the work of 
Stroud and Secrest (1966), are employed to automatically select the proper order of the 
quadrature rule. Additionally, to improve accuracy in a cost-effective manner, a graded 
subdivision of the element is incorporated, especially when ( is nearby. 
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The integration over the surface of the hole and fiber elements must be carried out 
using numerical integration. The complexity of the transient kernel function prohibits an 
analytic integration about the circumference of the fiber in a manner similar to the steady- 
state case. Therefore, numerical integration must be performed in both the longitudinal 
and the circumferential directions. An efficient integration scheme is adopted in which one- 
dimensional Gaussian integration formulas are applied independently in the two directions. 
This allows the subsegmentation and mapping, as well as the order of the Gaussian formulas 
to vary independently in the two directions so that the accuracy and the cost of the 
integration can be optimized. 

Turning next to category (2), one finds that again Gaussian quadrature is applicable, 
however, a somewhat modified scheme must be utilized to evaluate the weakly singular 
integrals. This is accomplished through element subsegmentation about the singular point 
so that the product of shape function, Jacobian and kernel remains well behaved. 

Unfortunately, the remaining strongly singular integrals of category (3) exist only in 
the Cauchy principal value sense and cannot, in general, be evaluated numerically, with 
sufficient precision. It should be noted that this apparent stumbling block is limited to the 
strongly singular portions, ss Fij and 33 Fge, of the kernel. The remainder of F*p, including 
tr Flj and tr F$ e , can be computed using the procedures outlined for category (2). However, 
as will be discussed in the next subsection, even category (3) 33 Fij and 33 Fee kernels can be 
accurately determined by employing an indirect ‘rigid body’ method originally developed 
by Cruse (1974). 

5.3.4 Assembly of Equations 

The complete discretization of the boundary integral equation, in both time and space, 
has been described, along with the techniques required for numerical integration of the 
kernels. Now, a system of algebraic equations can be developed to permit the approximate 
solution of the original quasistatic problem. This is accomplished by systematically writing 
(5.6) at each node on the outer surface, on the surface of the hole, and on the surface of 
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the fiber. The ensuing nodal collocation process, then produces a global set of equations 
given here in matrix form, for time N. 

For the nodes on matrix: 

[Gi,]{<£} - [Fjvf]{tio} 4- [G^K*//} - (5.7a) 

For the nodes on the fibers: 

[G 1 F ]{t%}-[F 1 F ]{u , *} = {D”} (5.76) 

where {D N } and {Dp } are the values of the convolution from the previous N - 1 time steps 
defined as 

~{D"} = £ {[G&+ l -"]{<3} - [< +1 - n ]{u2»} + [< +1 -"]{<h} - (5.8a) 

rj — 1 

and 

-{D?} = E {[GF +1 - n m} - [F£ +1 - n ){u n F }} (5.86) 

n=l 

In these equations, subscripts O, H and F are used to denote the quantities associated with 
the outer surface of the matrix, the surface of the hole containing the fiber, and surface of 
the fiber, respectively. Furthermore, the C Q/ g(£) term of each equation is lumped with the 
respective coefficient in the F 1 matrix. 

To calculate the coefficients of the singular points by the ‘rigid body’ technique, con- 
sider now, the first step. Thus, for N = 1, equation(5.7) becomes 

[G^lfro} ~ ^mH u o} + _ [Dm]{«h} = {0} (5.9a) 

[Gpjfrp} — [Dp]{up} = {0}. (5.96) 

At this point, the coefficients of the F matrices corresponding to the singular node of the 
equations has not been completely determined due to the strongly singular nature of the 
kernel function. To determine the values of these coefficients, we first decompose the Fij 
kernel into transient and steady-state parts. Following Cruse (1974) and, later Banerjee 
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et al (1986) in elastodynamics, the steady-state part of the coefficients can be calculated 
indirectly by imposing a uniform ‘rigid body’ generalized displacement field on the same 
body under steady-state conditions. The steady-state part of the singular coefficients 
are simply the summation of the non-singular, steady-state coefficients. The remaining 
transient portion of the coefficients axe non- singular, and hence can be evaluated to any 
desired precision. 

For a fiber that is perfectly bonded to the matrix, the displacement of the surface on 
the hole containing the fiber is equal to the displacement of the surface of the fiber, and 
the surface tractions are equal in magnitude but opposite in direction. 

{“&} = {“?} (5.10a) 

{<£} = -{*$} (5.10 b) 

Substitution of equation (5.10) into equations (5.7b) and (5.8b) and rearranging yields 

<<K} = -[ffM<«i5) -[»]{»?} (5.11a) 

{£>?) = Y. {[G? +, -"](lM + [F£ ,+, -'‘]<<*)) . (5.115) 

1 

where 

[H] = [ G ],]- 1 

Note, since the boundary integral equation for each fiber is independent of the equations 
for the other fibers, the inversion of [G^j can be reduced to a number of smaller inversions, 
one corresponding to each fiber. 

Equation (5.11a) can now be back substituted in equation (5.7a) to yield 

[<&]{#} - [Fh]{u%] - [F*]{ U £} = {D n } + [G}M{Z#} (5.12) 

where 

[F X ] = [G 1 h ][H}[FM + [F 1 h } 
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In a well-posed problem, at time A*, the set of global generalized nodal displacements 
and tractions will contain exactly (4 x P) unknown components (and P unknowns for heat 
conduction analysis). 

Then, as the final stage in the assembly process, equation (5.12) can be arranged to 
form 

[yl 1 ]{x Ar } = [5 1 ]{j/ v } + {^} + [Gy[fl']{D^} (5.13) 

in which 

{x N } unknown components of {u N } and {t^} 

{y^} known components of {/} and {t N } 

[A 1 ,^!? 1 ] associated matrices 

Note, the entire matrix [F x ] becomes part of [.A 1 ] since all the (generalized) displacements 
on the interface are unknown quantities. 

5.3.5 Solution 

A solution of equation (5.13) may be achieved for any time using a time marching 
algorithm. For the initial time N = 1, {Z? 1 } = {Z^} = {0} and equation (5.13) can be 
rewritten as 

[/l 1 ]{x 1 } = [5 1 ]{y 1 } (5.14) 

To obtain a solution for the unknown nodal quantities of this equation, a decomposition 
of matrix [A 1 ] is required. In general, [A 1 ] is a densely populated, unsymmetric matrix. 
The out-of-core solver, utilized here, was developed originally for elastostatics from the 
LINPACK software package (Dongarra et al, 1979) and operates on a submatrix level. 
Within each submatrix, Gaussian elimination with single pivoting reduces the block to 
upper triangular form. The final decomposed form of [A 1 ] is stored in a direct-access file 
for reuse in subsequent time steps. Backsubstitution then completes the determination of 
{z 1 }. Additional information on this solver is available in Banerjee et al (1985). 

After returning from the solver routines, the entire nodal response vectors, {u 1 } and 
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{t 1 }, at time At are known. For solutions at later times, a simple marching algorithm is 
employed. Assuming that the same set of nodal components axe unknown as in (5.14) for 
the first time step, equation (5.13) is reformulated as 

[A l ]{x 2 } = [5 x ]{y 2 } + {G 2 } + (5.15) 

Since the right-hand side contains only known quantities, (5.15) can be solved for {i 2 }. 

However, the decomposed form of [A 1 ] already exists on a direct-access file, so only the 

relatively inexpensive backsubstitution phase is required for the solution. 

The generalization of (5.15) to any time step N is simply equation (5.13) 

+ P"} + (5-16) 

in which the vectors {D N } and {£>£} contain summations which represent the effect of past 
events. By systematically storing all of the matrices and nodal response vectors computed 
during the marching process, surprisingly little computing time is required at each new 
time step. In fact, for any time step beyond the first, the only major computational task 
is the integration needed for form [G^] and [F N ]. Even this process is somewhat simplified, 
since now the kernels are non-singular. Also, as time marches on, the effect of events that 
occurred during the first time step diminishes. Consequently, the terms containing [G^] 
and [F n ] will eventually become insignificant compared to those associated with recent 
events. Once that point is reached, further integration is unnecessary, and a significant 
reduction in the computing effort per time step can be achieved. 

It should be emphasized that the entire boundary element method developed, in this 
section, has involved surface quantities exclusively. A complete solution to the well-posed 
linear uncoupled quasistatic problem with composite fibers can be obtained in terms of 
the nodal response vectors, without the need for any volume discretization. In many 
practical situations, however, additional information, such as, the temperature at interior 
locations or the stress at points on the boundary, is required. The next section discusses 
the calculations of these quantities. 
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5.4 Interior Quantities 

Once equation (5.16) is solved, at any time step, the complete set of nodal displace- 
ments and tractions are known. Subsequently, the response at points within the body can 
be calculated in a straightforward manner. For any point £ in the interior, the generalized 
displacement can be determined from (5.6a) or (5.6b) with Cp a = 6p a . Now, all the nodal 
variables on the right-hand side are known, and, as long as £ is not on the boundary nor 
on the interface between the fiber and the matrix. The kernel functions in (5.6) remain 
non-singular. However, when f is on the boundary or on the interface, the strong singular- 
ity in 8S Fp a prohibits accurate evaluation of the generalized displacement via (5.6), and an 
alternate approach is required. The apparent dilemma is easily resolved by recalling that 
the variation of surface quantities is completely defined by the elemental shape functions. 
Thus, for points on the outer boundary of the matrix, the desired relationship is simply 

«« «) = ( 5 - 17 « ) 

Where L w are the shape functions for the appropriate element and ( are the intrinsic coor- 
dinates corresponding to £ within that element. Obviously, from (5.17), neither integration 
nor the explicit contribution of past events are needed to evaluate generalized boundary 
displacements. 

In many problems, additional quantities, such as heat flux and stress, are also impor- 
tant. The boundary integral equation for heat flux, can be written 

*?(() = E { E [ / e% | 1 -"(X(C) - OMOrfs* WC))1 

n = 1 q= 1 

- E [/ W-"WO - OMOdS’WO)] 

9 _1 lvs» 

+ E [7 P + l - n (X(C) - ON„(OM^0)dSP(X{O) 

-til DpeV- n (X( 0 - OiV w «)M 7 (0)dSP(X(O)j ( 5 - 18 ) 
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where 


E2 ei (X(0-0 = ~k 


,8Gfa{X{ Q-0 


% 


D n 0ei (X(O-O = -k 


, 8Ffa(X( Q-0 


% 


(5.19a) 

(5.196) 


A similar equation can be written for the heat flux inside a fiber. This is valid for interior 
points, whereas, when £ is on the boundary, the shape functions can again be used. In this 
latter case, 


(5.20a) 

dL w ( C) ZJJ V_ 1 ..AVtN fc nnL\ 

—dT e “ ~'k ~dc qi (0 (5 - 206) 

which can be solved for boundary flux. Meanwhile, interior stresses can be evaluated from 


(5.21) 

(5.22a) 
(5.226) 

A similar equation could be derived if stress measurements are required inside a fiber. 
Equation (5.22) is, of course, developed from (5.6). Since strong kernel singularities appear 
when (5.22) is written for outer boundary points, an alternate procedure is needed to 
determine surface stress. This alternate scheme exploits the interrelationships between 
generalized displacement, traction, and stress and is the straightforward extension of the 


in which 




<«) = E { E f f E N+1 ~ n (x ( 0 - OMOdS’W 0) 

n=l < 2=1 *-^59 

- E [/ s< D%t'- -"(*«) -<)i.KW(0) 


l /3w 


E%! 1 - n (x ( 0 - O^(C)M 7 (0)rfS^(x(C)) 


*£[/.. 

- E [/, '■"(*«) -j.,} 


%(^(0-0 = 


—*«—+« I ^- + -55- 
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technique typically used in elastostatic implementation (Dargush and Banerjee, 1990). 
Specifically, the following can be obtained 


»j(0*S = LUOC 

(5.23 a) 

«i(0 + «ft(0) = -/Wtf £«(<)«& 

(5.236) 

dxj v"(e)= 9Lu ^K^ 

d( d( ^ 

(5.23c) 


in which is obviously the nodal temperatures, and, 

Dij M = \6ij6ki + 2nSik6ji. 

Equations (5.23) form an independent set that can be resolved numerically for <r£[(£) and 
u ij(0 completely in terms of known nodal quantities and without the need for kernel 
integration nor convolution. Notice, however, that shape function derivatives appear in 
(5.23c), thus constraining the representation of stress on the surface element to something 
less than full quadratic variation. The interior stress kernel functions, defined by (5.22) 
are also detailed in Appendix B. 
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6. NONLINEAR MATERIAL FORMULATION 

6.1 Introduction 

Due to the material discontinuity in a composite structure, large localized stress gra- 
dients develop on the fiber-matrix interface when a specimen is loaded. These large stress 
gradients induce plastic yielding in a composite about the fibers. It is therefore, impor- 
tant to include these material nonlinearities in an analysis of a composite structure loaded 
beyond its elastic limit. As a first step in this analysis, a formulation to account for the 
plastic yielding of the composite matrix is developed. The nonlinear effects are included 
in the boundary element formulation through a volume integral. 

In this chapter, new methods are developed for the efficient evaluation of the volume 
integrals which appear in the integral formulations. An elastoplastic-fracture constitutive 
law is used to model the nonlinear behavior of ceramic composite material, and an iterative 
control algorithm is implemented for the solution of the nonlinear system. 

6.2 Integral Equation Formulation for Elastoplasticity 

The governing equation of plasticity can be written in the following form [Banerjee 
and Butterfield, 1981]: 

(A + n)u j}ij + fiiiijj = <r c ijd (6.1) 

«»J = 1.2,3 

where A and n are Lame constants, u, is displacement rate and <7£ is an initial stress rate 
(or a corrective stress rate) resulting from the nonlinearities present in the plastic domain. 
For the purpose of the plasticity algorithm the initial stress rate is defined as 

** ij = Dijkl^kl ~ D% kl Skl 

where D? jkl and D^ kl is the elastic and elastoplastic constitutive tensors, respectively, and 
e k i is the strain rate. 
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The boundary integral equation for the displacement at a point £ in the matrix can be 
expressed as 

Ciiii)Mt)= [ [G%{z,Oi?(x)- F!?(x,Ou?(x))dS( X )+ f B^(z,0&? k (z)dV(z) 

J s J V 

+ £/ [G%(x,Oi”(x) - F’fix^ix^dS-ix) (6.2) 

n= 1 dS" 

where 

G^j, Fff, and B& k are the fundamental solutions of a ceramic matrix of infinite extent, 

Cij are constants determined by the geometry at f , 
iiiji are the displacement and traction rates, 

S, S n are the surfaces of the outer boundary and the nth hole (encompassing 
the fiber), respectively and 
N is the number of individual fibers. 

Superscript O and H identify the quantities of the outer surface of the matrix and the 
quantities on the surface of the hole, respectively. 

A similar equation can be written for the displacement rate of a point £ in the nth fiber 
as 

GfjUiiO- f [GF j (x,0i?{x)-F t F j {x,0^(*)]dS n (x) (6.3) 

JS n 

where 

Gfj, FF are the fundamental solutions of the nth fiber, 
are constants determined by the geometry at £, 
uf , if are displacement and traction rates associated with the nth fiber 
S n the surface of the nth fiber 
Note, each fiber may have different material properties. 

For a perfectly bonded fiber-matrix interface the displacement rates of the matrix and 
the fiber are equal and the traction rates are equal and opposite: 

uf(x) = uf(x) (6.4 a) 

<?(*) = -<f (*) ( 6 - 46 ) 
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As in the derivation of the integral equation for elastostatic fiber composite analysis of 
Chapter 2, the Poisson’s ratio of the fiber and the Poisson’s ratio of the matrix is assumed 
to be equal. This leads to the following relation: 

Fg(x,t) = -Fg(x,t) ( 6 . 5 ) 

Finally, using relations (6.4) and (6.5), equation (6.3) can be backsubstituted in equa- 
tions (6.2) to yield the integral equation for the displacement rate at a point £ in a fiber 
composite structure in which material nonlinearities axe present in the matrix: 

C<j(t) = J s [G%(x,Oi?(x) - F t y(x,Ou?(*)]dS(x) + ^B^(z,^? fc (*W*) 

N r 

+ £/ G ij (x,Oif(x)dS n (x) ( 6 . 6 ) 

n =i J s » 

where 

G ij (x,Z) = G? j (x,t)-Gf j (x,Z) 

Cij({) are constants dependent on the geometry for a point f on the outer boundary, 
and 

Cij(£) = Sij for a point £ in the interior of the composite. 

6.3 Cylinder Volume Cells 

The volume integral in equation (6.6) is evaluated assuming a variation based on a 30 
node shape function over a cylindrical volume cell (Figure 6.1). The initial stress rates have 
a variation through the volume cell which axe quadratic in the longitudinal direction, linear 
in the radial direction, and trigometric in the circumferential direction. The variation in the 
circumferential direction is based on a 5-noded circular shape function shown in Figure 6.2. 
It is similar to the 3-noded circular shape function used for the displacement and traction 
rates on the fiber/hole, however, in addition to the constant, sine and cosine terms, the 
5-noded function contains sin 29 and cos 2 6 terms. Hence, the shape function can model 
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initial stress rates which have two peak values in the circumferential direction. In two- 
dimensional numerical experimentation, this was observed as a minimum requirement for 
reasonable approximation of the actual initial stress rate variation. 

The integration is carried out analytically about the circumference of the volume cell 
to keep numerical costs down. However, due to the complexity of the resulting integrand, 
numerical integration was employed in the radial and longitudinal directions. 

6.4 Interior Stress Rate Equations 

The solution algorithm in a plasticity analysis requires the determination of stress (and 
strain) rates at nodes of the volume cell. The integral equation for stress (and strain) rates 
in the composite matrix are derived from the displacement rate integral equations (6.2) 
through the application of the strain-displacement equation and the nonlinear stress-strain 
relations, i.e., 

ctij = Df jkl e kl - <r% ( 6 . 7 ) 

Nonlinear effects are incorporated in the stress rate integral equation via a volume integral. 
This volume integral is strongly singular, however, and special considerations must be given 
to this integral as described in Henry and Banerjee (1987). 

The integral equation for stress rates, however, are strongly singular for a point on 
the fiber-matrix interface. In this instance the stress rates are best determined using a 
boundary stress calculation similar to equation (2.13). However, in a nonlinear analysis the 
nonlinear stress-strain rate constitutive relationship involving initial stress rates [Henry, 
1987] must be used. The calculation is very efficient because the calculation does not 
require surface nor volume integration. 

6.5 Plasticity-Fracture Constitutive Model for Ceramic Composites 

An elastic-plastic strain hardening-fracture model (Chen, 1975) is used to model the 
inelastic behavior of the ceramic material. The yield function describing the inelastic 
behavior is bounded by two surfaces: the elastic limit or initial discontinuous surface; and 
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a failure surface. When the stress level violates the elastic limit, plastic yielding occurs 
with strain hardening. At this point unloading would behave linear-elastically, however, 
permanent deformation will have occurred. A new elastic limit surface, uniquely defined by 
the invariants p and 3 2 (described below), corresponds to the highest stress state achieved 
in prior loading history. Loading can continue until the failure surface is reached. At this 
time the material will fracture and the stress will be reduced to zero at this point. 

An important feature that this model exhibits to accurately represent ceramic material 
behavior is its ability to sustain stress in compression many times greater than in tension. 
Further, this phenomenon is magnified during plastic hardening since both the isotropic 
and kinematic hardening rules axe employed. The strain-hardening model is illustrated in 
Figures 6.3 and 6.4. 

The loading function shown in Figure 6.4 has the following form: 
in the compression domain: 


f( n \ h + * Ap r 
~ 1 -\Bp~ T 


in the tension and tension-compression domains: 


f(a ., _ Ji-y 2 +k*p _ T 2 

H tj) ~ \-\Bp 


where 

p = ((7 X + <T y + <T Z )/ 3, 

h - I [(<7® - Oy) 2 + (<T y ~ z) 2 + {? z ~ <^x) 2 } + *ly + + a lzl 

A - a ° t 7 u~ a ? t ? 


(6.8a) 


(6.86) 


7 u '0 . 


r is a function of hardening, and 

Ao, t 0 , A u , r u are material constants which can be obtained from uniaxial and biaxial 
tests. 
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By employing the normality condition the final form of the constitutive equation relating 
the stress increment to the strain increment can be derived. 

df df 


Okl = 


Dim ~ PD 


r>e 

klmn d(T bn W 
(J °mn uo qr 


( 6 . 9 ) 


where 


* = (%-$-)* + &d, 


w'aui.i 


Di Uj i s the elastic constitutive matrix relation 

H' is the current slope of the equivalent stress-equivalent plastic strain hardening curve, 
and 

o e = y/J is the definition of equivalent stress 


6.6 Assembly of Equations 

The discretization and integration of the integral equations for displacement and stress 
rates follows the procedure outlined in Chapter 2. After the integration of the kernel 
functions, the next critical step in the formulation is the assembly of the system equations. 
Here, efficiency is of utmost importance. The approach to writing an efficient algorithm 
is to keep the number of system equations to a minimum by eliminating all unnecessary 
unknowns from the system. The strategy is to retain in the system only traction variables 
on the fiber-matrix interface. This is in contrast to a general multi-region problem where 
both displacement and tractions axe retained on an interface. The elimination of the 
displacement on the interface is achieved through a backsubstitution of the fiber equations 
in the system equations which are made up exclusively from equations written for the 
composite matrix (on the outer surface and on the surface of the holes). The procedure is 
described below. 

Equation (6.6) is used to generate a system of equations for nodes on the outer surface 
of the composite matrix and for nodes on the surface of the holes containing the fibers. 
Written in matrix form we have 


On the Matrix Outer Surface: G M t° - F M u° + Gt H + B M <7 C = 0 (6.10a) 
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On the Matrix Hole Surface: G M t° — F M u° + Gt H + B M o- c = Iu H (6.106) 

where 

t° and u° are traction and displacement rates 
t H and u H are traction and displacement rates 
<r c is the initial stress rate vector 
I is the identity matrix 

G m and F M matrices contain coefficients from the integration over the outer boundary. 
G matrix contains coefficients integrated about the fiber /hole 

B m matrix contains coefficients from the integration over the volume cells 

Our goal is to eliminate u H from the system. To this end, equation (6.3) is written for 
every node on a fiber, collocating slightly outside the boundary of the fiber [at a distance 
of (1.25)*(fiber radius)] where c£(0 = 0. 

F F2 u F = G F2 t F 

Superscript F2 identifies the equations written at points located slightly outside the bound- 
ary of the fibers. 

Noting u M = u F and t H = -t F we have 

F F2 u h = — G F2 t H (6.11) 

Post multiplying equation (6.10b) by the F F2 matrix in equation (6.12) yields 

pF2 G M jO _ F F2pM^O + pF2QjH + pF2gM ( j.C _ pF2^H (6.12) 

Equation (6.11 ) can now be set equal to equation (6.12) and the final form of the system 
is derived. 

On Outer Surface: G M t° - F M u° + Gt H + B m <t c = O 

On Hole: F F2 G M t° - f F2 F m u° + (F F2 G + G F2 )t H + F F2 B m <t c = O (6.13) 
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At every point on the outer surface, either the traction or the displacement rate is 
specified and on the surface of the hole only the traction rates axe retained. Therefore, 
the number of displacement rate equations in the system are equal to the final number 
of boundary unknowns. The initial stress rate vector, however, is unknown and must be 
determined by a quasistatic, iterative solution algorithm. This requires the use of stress 
rate equations and a constitutive law. The stress rate equation can be written in matrix 
form as 

<r = G a i° - F'u 0 + G't 11 - F a u H -h B a & c (6.14) 

Note, the displacement rate vector u H for the hole is known from equation (6.10b) 
and could be eliminated from equation (6.14). For numerical efficiency, however, a two 
step procedure is used in which u H is first determined using equation (6.10b) and then 
substituted in equation (6.14). Nevertheless, for illustration purposes, a backsubstitution 
of equation (6.10b) into equation (6.14) will be implied. Hence, the displacement rate and 
stress rate equations (6.13) and (6.14) are assembled by collecting the known and unknown 
values of traction and displacement rates and their coefficients together. The final system 
equations can be cast as: 

A b x = B b y + C b <r c (6.15a) 

<r = A^x + B^y+C* 7 ^ (6.156) 

where x is the vector of unknown variables at boundary and interface nodes; y is the vector 
of known variables; <r c is the vector of initial stress rates; A b , B b , C b are the coefficient 
matrices of the boundary (displacement rate) system equations (6.13); and A a ,B a ,c a axe 

the coefficient matrices of the stress rate equation (6.14). It should be noted that A b is a 

> 

square matrix. Furthermore, in a substructured system the matrices A b and B b are block 
banded while matrices C^A", B a , and C° are block diagonal. 
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6.7 Nonlinear Iterative Solution Algorithm 

The algorithm described below provides the solution of the system equations given by 
(6.15). This requires complete knowledge of the initial stress rate distribution b c within the 
yielded region that is induced by the imposition of the boundary loading. Unfortunately 
the nonlinear initial stress rates are not known a priori for a particular load increment and 
therefore an iterative process must be employed within each loading stage. 

An important feature incorporated in the iterative algorithm of the present work is 
an acceleration scheme based on the initial stress rates generated by the past history. In 
this procedure, the path followed by the previous load increment is used to extrapolate 
the initial stress rates at the beginning of the current increment before the start of the 
iterative operations. This results in substantial reduction in computer time. This proce- 
dure is graphically illustrated in Figure 6.5 for a simple tension problem with a variable 
hardening parameter. The initial loading from point A generates initial stress rates which 
are distributed during the iterations to arrive at the solution at B. Upon loading from 
state B the extrapolated path BC is followed. At C the resulting initial stress rates are 
distributed to reach state D. 

The steps of the incremental plasticity algorithm is described below: 

1. Obtain the elastic solution x for an arbitrary increment of boundary loading y from 

A b x = B b y (6.16a) 

and determine the stress rates at the nodes of the volume cells from 

b = A. a i + B°y (6.166) 

2. Scale the elastic solution such that the node with the highest stress rate is at yield. 
In the case of nonproportional loading, this onset of yielding cannot be reached by 
scaling. Instead, incremental loads must be applied until yielding is reached. 

3. Impose a small load increment y, (usually less than 5% of the yield load) plus the 
estimated value of the initial stress rates accumulated from the previous load step and 
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evaluate x and a using 


A b x = [B b y + C b «r c ] (6.17a) 

b = A a x + B°y + C°b c (6.176) 

where a c is the estimated initial stress rates. If no prior plastic history exists the value 
of the estimated initial stress rates axe zero. 

4. Accumulate all incremental quantities of stress, traction, and displacement rates and 
use Eq. (6.7) to calculate the strain for this increment. 

5. Evaluate the current constitutive matrix using the new stress rate history and calculate 
the current initial stress rates via 

&fj = (Df jkl - D% kl )£ kl (6.18) 

Accumulate the initial stress rate history of this load increment for use as the estimated 
initial stress rate of the next load increment in step 3. 

6. If the current increment of initial stress rates computed in Eq. (6.18) is greater than 
a prescribed tolerance (normally 0.005 times the yield stress) then calculate the incre- 
mental quantities (x and a) due to these rates using 

A b x = C b (7 c and b = A a x 4- C°b c 

and return to step 4 for the next iteration. If the value is less than the prescribed 
tolerance, go to step 7. Note the boundary loading is zero for this calculation. If the 
number of iterations is greater than a specified limit (usually 50) the system is assumed 
to have reached the state of failure. 

7. Return to step 3 and apply the next load increment and the accumulated initial stress 
rates from this load step. (If the size of the load increment changes the estimated 
initial stress rates should be scaled proportionally.) 

Any residual initial stress at the end of the iteration is carried forward and applied to 
the system with the next load increment. 
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%(()»<«) . /;<*(..*?(«) - C(..!).?( t )l J S W + 1 

N r 

+ E y 5 . «?£(*. ow - 06, w (*)irfs" ( ,) 



c£ui(0 = / lc5(x,0‘T(*)- F y(*>0“r(*)] dSB (*) 

J ./s* 



Cy(0 = jf - f^ (l ^) u p( r )]dS(r) + jf^OiSlW) 

+ E/ 6 y (*. ()*'«"<*)«"(*) 

n = l ‘ /5 ’' 

where 

6 y (c,0»cg(*.0-c y («.0 

»',i = 1,2,3 


Fig. 6.1 Integral Equation Formulation for the Nonlinear Analysis of Fiber Composite 
Materials. A 30-Node, Cylindrical Volume Cell is Generated Automatically. 
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CIRCULRR SHRPE FUNCTION VRLUE 
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Fig. 6.2 Value of the 5-Noded Circular Shape Function about the Circumference of 
the Cylindrical Volume Cell 



Fig. 6.3 Typical Stress-Strain Curve in Equivalent Stress-Strain Space 
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Fig. 6.4 Loading Function in Biaxial Principal Stress Space 



Fig. 6.5 Accelerated Iterative Plasticity Algorithm 
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7. COMPUTER PROGRAM DEVELOPMENT 

7.1 Introduction 

The goal of the ‘BEST-CMS’ computer program developed for ceramic composites is 
the accurate and efficient implementation of the formulation described in Chapters 2-6. 
Of equal importance is the degree of generality required in the definition of component 
geometry, loading and material properties. This is necessary if the program is to be 
applicable to real problems in the aerospace industry. 

For this reason the ceramic composite formulation has been implemented in the three- 
dimensional boundary element computer code ‘BEST3D’ (Boundary Element Stress Tech- 
nology - Three-dimensional) which was developed for NASA by Pratt and Whitney and 
SUNY/Buffalo under contract NAS 3-23697. Since its development, BEST3D has been 
proven to be a highly accurate and numerically efficient boundary element program. 

The development of the computer program ‘BEST-CMS’ is discussed in the following 
sections. 

7.2 Program Structure 

BEST-CMS is designed to be a fully general ceramic composite analysis system em- 
ploying the boundary element method. The program is written using standard FORTRAN 
77. Development has been carried out at SUNY/Buffalo on an HP 9000 minicomputer, a 
SUN-4 workstation, and a SUN Sparcstation-1. The nature of the method is such that, 
for any realistic problem, not all required data can reside simultaneously in core. For this 
reason extensive use is made of both sequential and direct access scratch files. 

The program first executes an input segment. After the input has been processed, 
the surface integrals are calculated and assembled into the set of system equations using 
specified boundary conditions, followed by the fiber assembly and the inclusion of the fiber 
equations in the general system. The system matrix is then decomposed and saved on 
disk, followed by the calculation of the solution vector. The full displacement and traction 
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(and/or temperature and flux) solution on each boundary element and fiber element is 
then reconstructed from the solution vector. In a time dependent problem the process of 
constructing the load vector for the system equations is repeated at each time step, but 
the integration, formation and decomposition of the system matrix are done only once. 
However, in an analysis with nonlinear fiber-matrix interface connections, the assembly 
and decomposition must be carried out several times. 

The current program limits are as follows: 

- 20 time points (solutions for different loadings) 

- 15 generic modeling regions (GMR) 

- 600 elements (300 elements per GMR) including fiber elements 

- 100 fiber elements per generic modeling region 

- 500 fiber elements per problem 

- 2500 modeling nodes 

- 1200 source points (600 source points per GMR) 

Various aspects of the computer program are discussed below. 

7.3 Program Input 

The input for BEST-CMS is free field. Meaningful keywords are used to identify data 
types and to name particular data sets. The input is divided into five types: 

1. Case Control Cards 

The case control cards define global characteristics of the problem. In addition to the 
problem title, the times for multiple time steps are defined. The reading or writing of 
restart data is also defined at this point. The restart facility allows one to change 
the arrangements to fibers without recalculating the various coefficients. 


2. Material Property Definition 

The material property input allows the definition of material properties for a variety 
of materials. The Young’s modulus can be prescribed in tabular form for a user-defined 
set of temperatures. Temperature independent values of Poisson’s ratio are also defined. 

3. Geometry Input 

Geometry input is defined one GMR (generic modeling region, or subregion) at a time. 
To initiate the input, a tag is provided to identify the GMR, a material name and reference 
temperature are defined to allow initialization of material properties. 

The next block of geometry input consists of the Cartesian coordinates of the user 
input points for the outer surface geometry definition of the composite matrix, together 
with identifiers (normally positive integers) for these geometric nodes. 

Following the definition of an initial set of nodal points, the surface connectivity of the 
outer surface of the composite matrix is defined through the input of one or more named 
surfaces. Each surface is made up of a number of elements, with each element defined 
in terms of several geometric nodes. Three sided elements, defined using six rather than 
eight geometric nodes, are used for mesh transition purposes. The terms quadrilateral and 
triangle are normally used to refer to the eight and six noded elements, although the real 
geometry represented is, in general, a nonplanar surface patch. Nine noded elements are 
made available by adding a central node to the eight noded elements. 

Over each element the variation of displacement and traction (and/or temperature and 
flux) can be defined using either the linear, quadratic or quartic shape functions. Linear 
and quadratic elements can share a common side, which is then constrained to have linear 
displacement and traction (and/or temperature and flux) variation. Quadradic/quartic 
mi xtures assumes a quadratic variation, and linear /quartic mixtures assumes a linear vari- 
ation. 

Finally an option is available to allow quadratic and quartic functional variation to be 
used in conjunction with linear geometry (4 or 3 nodes). In this case the program generates 
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the additional nodes automatically at mid-point of the sides. 


various element types are summarized below. 

Geometry 


Surface Element Type Nodes 

Linear Quadrilateral 3, 8 or 9 

Linear Triangle 4, 6 or 7 

Quadratic Quadrilateral 4, 8 or 9 
Quadratic Triangle 3 or 6 

Quartic Triangle 3, 8 or 9 

Quartic Quadrilateral 3 or 6 


The characteristics of the 


Displacement/Traction Nodes 
(and/or Temperature/Flux Nodes) 

4 

3 

8 or 9 

6 

13 or 15 
17 or 25 


Following the definition for the composite matrix outer surface, the embedded fibers 
are then defined. These are defined as curvilinear line elements with a prescribed radius 
of the cross-section. The fibers are generally straight, however as noted, curved fibers are 
also allowed. The user first defines the nodal coordinates of the centerline of the fiber. 
Thereafter, the radius and the fiber connectivity is defined. Linear and quadratic elements 
are available for both geometry and functional variation, however, quadratic functional 
variation over linear geometry is not presently available. The various options for the fiber 
elements are summarized below. 

Geometry Displacement/Traction Nodes 
Fiber Element Type Nodes (and/or Temperature/Flux Nodes) 


Linear- Linear 
Quadratic- Linear 
Quadratic- Quadratic 


2 

3 

3 


2 

2 

3 


Note only the surface of the fiber needs to be defined, i.e., the hole in the composite matrix 
which encompasses the fiber does not have to be explicitly defined. 

4. Interface Conditions 

The interface input describes the connection of surfaces or elements of one region 
to another and between the matrix and the fibers. Special types of fiber-matrix interface 
conditions which are available presently include fully-bonded and sliding contact (including 
coulomb friction) and springs. 
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5. Boundary Condition Input 

The final input section provides for the definition of boundary conditions, as functions 
of both position and time. Data can be input for an entire surface, or for a subset (ele- 
ments or nodes) of a surface. Input can be in global coordinates, or can define rollers or 
pressure (or flux) in the local coordinate system. Input simplifications are available for 
the frequently occurring cases of boundary data which is constant with respect to space 
and/or time variation. Each boundary condition set can be defined at a different set of 
times. 

7.4 Surface Integral Calculation 

Following the processing of the input data, the surface integrals are evaluated numeri- 
cally. This is the most time consuming portion of the analyses. In BEST-CMS the results 
of these integrations are stored as they are calculated, rather than being assembled into 
the final equation system immediately. Although this is somewhat more costly in terms of 
storage and CPU (central processing unit) time, it has led to much greater clarity in the 
writing of BEST-CMS. In addition, it provides much greater flexibility in the implemen- 
tation of various restart and boundary condition options. 

The calculations proceed first by GMR (generic modeling region), then by source point 
(the equation being constructed) and finally by surface element and fiber element. The 
results for each source point element pair are written to disk. All of the calculations are 
carried out and stored in the global (Cartesian) coordinate system. 

The integration of the BEM equations is the most complex part of the code. In this 
process either singular or nonsingular integrals can be encountered. The integrals are 
singular if the source point for the equations being constructed lies on the element being 
integrated. Otherwise, the integrals are nonsingular, although numerical evaluation is still 
difficult if the source point and the element being integrated are close together. 

In both the singular and nonsingular cases Gaussian integration is used. The basic 
technique is developed in Banerjee and Butterfield, 1981. In the nonsingular case an 
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approximate error estimate for the integral was developed based on the work of Stroud 
and Secrest (1966). This allows the determination of element subdivisions and orders of 
Gaussian integration which will retain a consistent level of error throughout the structure. 
Numerical tests have shown that the use of 3, 4 and 5 point Gauss rules provide the best 
combination of accuracy and efficiency. In the present code the 4 point rule is used for 
nonsingulax integration, and error is controlled through element subdivision. The origin of 
the element subdivision is taken to be the closest point to the source point on the element 
being integrated. 

If the source point is very close to the element being integrated, the use of a uniform 
subdivision of the element can lead to excessive computing time. This frequently happens 
in the case of aerospace structures, due either to mesh transitions or to the analysis of 
thin walled structures. In order to improve efficiency, while retaining accuracy, a graded 
element subdivision was employed. Based on one-dimensional tests, it was found that 
the subelement divisions could be allowed to grow geometrically away from the origin of 
the element subdivision. Numerical tests on a complex three-dimensional problem have 
shown that a mesh expansion factor as high as 4.0 can be employed without significant 
degradation of accuracy. 

In each case of singular integration (source point on the elements being integrated) 
the element is first divided into subelements. The integration over each subelement is 
carried out using a Jacobian transformation in mapping. This coordinate transformation 
produces nonsingulax behavior in all except one of the required integrals. Normal Gauss 
rules can then be employed. The remaining integral (that of the traction kernel FiJ times 
the isoparametric shape function which is 1.0 at the source point) is still singular, and 
cannot be numerically evaluated with reasonable efficiency and accuracy. Its calculation is 
carried out indirectly, using the fact that the stresses due to a rigid body translation axe 
zero (Lachat and Watson, 1976). It has been found that subdivision in the circumferential 
direction of a two-dimensional surface element is required to preserve accuracy in the 
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singular integration of the outer surface. A maximum included angle of 15 degrees is used. 
Subdivision in the radial direction has not been required. 

The integrals required for calculation of displacement, stress, temperature, and flux at 
interior points are of the same type as those involved in the generation of the system equa- 
tions, except that only nonsingular integrals are involved. If the source point involved is 
located on the surface of the body, then numerical integration is not required. Instead, the 
required quantities are calculated using the displacements and tractions (or temperature 
and flux) on the element (or elements) containing the source point, as discussed in Section 
2.5. 

7.5 System Matrix Assembly 

The first step in the assembly process is the reduction of the rectangular matrix of 
F integrals to a square matrix. This matrix is the prototype of the system matrix. The 
columns of the matrix are transformed or replaced, as required by the boundary conditions, 
as the assembly process proceeds. 

The next step in the process is the incorporation of the fiber equations in the system. 
As was described in detail for each specific analysis in Chapters 2-6, the fiber assembly 
consists of a fiber by fiber matrix multiplication and backsubstitution. The backsubstitu- 
tion minimizes the number of equations required in the system by eliminating some of the 
unk n own quantities on the fibers. 

A key problem in the entire process is the proper definition of appropriate coordinate 
systems, on a nodal basis. This is a problem common in any direct boundary element 
method which treats structures with nonsmooth surfaces. It arises because the tractions 
at a point axe not uniquely determined unless the normal direction to the surface varies 
continuously at the point in question. 

The original surface integral calculations are all done in global coordinates. If the 
displacement (or temperature) boundary condition is specified at a given node, in global 
coordinates, then no new coordinate system definition is required. It is only necessary to 


72 


keep track of the subset of elements, containing the given node, on which the fixed dis- 
placement (or temperature) is to be reacted. However, if a displacement (or temperature) 
is specified in a nonglobal direction at a given node, then a new nodal coordinate system 
must be defined and, potentially, updated as further boundary conditions axe processed. 
The associated nonzero reactions must then be expressed in the new coordinate system. 

Following this preparatory work, the final assembly of the system equations is carried 
out. It is performed in three major steps: 

1. Transformation of the columns of the matrices to appropriate local coordinate systems 
and incorporation of any boundary conditions involving springs. 

2. Incorporation of compatibility and equilibrium conditions on interfaces. 

3. Application of specified displacements and tractions (and/or temperatures and fluxes). 

Two particular features of the equation assembly deserve special comment. First, 
in multi-GMR problems the system matrix is not full. Rather it can be thought of as 
consisting of an NxN array of submatrices, each of which is either fully populated or 
completely zero. Only the nonzero portions of the system equations are preserved during 
system matrix assembly. In order to improve the numerical conditioning of the system 
matrix for the solution process, the columns are reordered so that the variables from the 
two regions, lying on the same interface, are as close together as possible. 

Second, rather than simply assembling an explicit load vector at each time point in the 
solution process, load vector coefficient matrices axe assembled and stored. These allow 
the updating of the load vector at any required time point simply by interpolating the 
time dependent boundary conditions and performing a matrix multiplication. A Similar 
process is used in the calculation of interior and boundary stresses. 
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7.6 System Equation Solution 

The solver employed in BEST-CMS operates at the submatrix level, using software 
from the LINPACK package (Dongarra, 1979) to carry out all operations on submatrices. 
The system matrix is stored, by submatrices, on a direct access file. The decomposition 
process is a Gaussian reduction to upper triangular (submatrix) form. The row operations 
required during the decomposition are stored in the space originally occupied by the lower 
triangle of the system matrix. Pivoting of rows within diagonal submatrices is permitted. 

The calculation of the solution vector is carried out by a separate subroutine, using the 
decomposed form of the system matrix from the direct access file. The process of repeated 
solution, required for problems with multi-time steps, is highly efficient. 

7.7 Nonlinear Solution Process 

The nonlinear solution algorithm starts with an elastic analysis (with linear elastic 
spring interface or bonded connection between the matrix and the fiber) of the problem for 
the first loading increment (complete with the specified boundary and body force loading). 
At the end of the elastic increment the state variables are calculated and the nonlinear 
constitutive relations are established. The difference between the actual values and the 
elastic (or previous) values is estimated. A new equation system is assembled with the 
calculated nonlinear constitutive relations. The process is essentially repeated until the 
constitutive equations yield values that are negligibly different from the previous step. A 
new loading increment is taken and the process is repeated for each subsequent increment. 

7.8 Output Description 

The output from BEST-CMS is relatively straightforward. It consists of ten sections, 
as follows: 

1. Complete echo of the input data set. 

2. Summary of case control and material property input. 

3. Complete definition for each GMR, including all surface fiber nodes, surface and fiber 
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elements. 

4. Complete summary for each interface and boundary condition set, including the ele- 
ments and nodes affected. 

5. Boundary solution (on an element basis), including displacements and tractions (and/ 
or temperature and fluxes) at each node of each element. 

6. The resultant load on each element and on the entire GMR is calculated and printed. 

7. Solution for the displacements and tractions (and/or temperatures and fluxes) at the 
Fiber-Matrix composite interface (on an element basis). 

8. Displacement, stress and strain (and/or temperature and fluxes) on a nodal basis, at 
all surface nodes, for each GMR. 

9. Displacements (and/or temperatures) at interior nodes. 

10. Stresses at interior nodes. 
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8. EXAMPLES OF FIBER COMPOSITE ANALYSIS 

8.1 Introduction 

In this section a number of examples are presented to verify and demonstrate the 
applications of the ceramic composite formulation of the previous chapters. 

In the mesh diagrams of the models containing the fibers, a double line is used to 
indicate the centerline of the fiber elements. The length of these elements are shown in 
proper proportion for the three-dimensional views, however, the radii of the fibers are not 
indicated on these diagrams. The double line is a symbolic representation of the fiber 
elements and does not in any way indicate the diameter of the fiber. Refer to the example 
description for the values of the radii. 

Throughout this section consistent units are used in the definition of the examples. 
This means all lengths are defined in the same units and the tractions and the elastic 
moduli are defined in terms of these lengths as force/length 2 . No confusion should arise 
since the results are reported as non-dimensional quantities. 

8.2 Cube with a Single Fiber 

The first test of the formulation is on a unit cube with a single fiber of radius 0.1 
through the center of the cube. The cube is subjected to tension and shear in the direction 
parallel and perpendicular to the fiber. The cube has a modulus of 100.0 and a Poisson 
ratio of 0.3. Consistent, units are used for all information described in this problem. A 
fiber with two different moduli of 1,000 and 10,000 is studied. The Poisson ratio of the 
fiber is assumed to be the same as that of the cube. 

The problem is analyzed by both the present formulation and by a full three-dimensional 
multi-region BEM approach. As shown in Fig. 8.2.1, the model for the fiber formulation 
consists of fourteen quadratic boundary elements and the fiber contains three quadratic 
fiber elements. The two-region, three-dimensional model shown in Fig. 8.2.2 contains 
twenty quadratic boundary elements in the first region and sixteen in the second. Note 
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9-noded elements axe used in describing the fiber and hole to accurately capture the curvi- 
linear geometry. 

In Fig. 8.2.3, the profile of the end displacement of the cube under a uniform normal 
traction of 100.0 (in parallel with the fiber) is shown. The present formulation is in good 
agreement with the full three-dimensional results for Ef/E = 10. For the case Ef/E = 100, 
the fiber formulation exhibits greater stiffness than the 3-D results. This difference is 
contributed to the way the load is distributed from the fiber to the composite matrix. In 
the full 3-D model, the applied traction and the resulting reactions at the fixed end act 
directly on the end of the fiber. In the composite formulation, the fiber is assumed not to 
intersect the boundary surface and therefore the fiber is moved back slightly from the end 
of the cube. The load is therefore transferred through the composite matrix to the end of 
the fiber and to its sides in a manner that is slightly different from the full 3-D analysis. 

In Fig. 8.2.4, the stress distribution through the center of the cube (from A to B 
as indicated in the figure) is shown. Again the results are very good for Ef/E = 10, and 
deviates slightly from the full 3-D results in the second case. 

In Figs. 8.2.5 and 8.2.6, the lateral displacements along the side of the cube axe shown 
for a cube subjected to a shear traction of 100. For the case of applied shear perpendicular 
to the fiber (Fig. 8.2.5), the results for both the fiber and full 3-D model show good 
agreement. Once again a slight deviation is observed for Ef/E = 100. In the case of 
the shear traction in the plane of the fiber (Fig. 8.2.6) the fiber has little effect on the 
displacement (as anticipated) and all results fall in close proximity. 

8.3 Lateral Behavior of a Cube with Multiple Fibers 

Existing methods of analysis of composite material based on mechanics of materials 
have been relatively successful in predicting the behavior of composite material for loading 
in the longitudinal direction. The properties perpendicular to the direction of the fibers 
are not so readily predictable by present means. The focus of the present example concerns 
this lateral behavior. 
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Four cubes (Fig. 8.3.1) with one, two, five and nine fibers are fixed with a roller 
boundary condition on one side and subjected to a uniform traction, perpendicular to the 
fibers. The material properties, given in consistent units, are 

E libsr = 1000 E matrix = 100 

l ,J iber = 0.3 /i matrix _ 0 3 

For the cube with one and two fibers, the boundary mesh consists of two quadratic 
surface elements on each lateral side and four elements on the top and bottom. For the 
cubes with five and nine fibers, one additional element was added to the side with the 
applied traction and to the side with the roller boundary condition. The top and bottom 
faces contain six elements to match the pattern of the sides. In all cases, each fiber 
contained three one-dimensional quadratic elements. 

The profile for the end displacement of the cube with one fiber and five fibers are 
shown in Figs. 8.3.2 and 8.3.3. The results are seen to be in good agreement with the 
two-dimensional results. The 2-D results are approximations since plane stress is assumed. 
The 3-D solutions for the one fiber are within 2% error of the 2-D solution and within 3% 
for the case of 5 fibers. 

Also shown in Fig. 8.3.4 are the average end displacements for the one, two, five and 
nine fibers. Results show good agreement with 2-D results. For one, two and five fibers, 
the solutions are within 2% error of the 2-D results and 6% for the case of nine fibers where 
the volume ratio of the fibers to the total volume is 28.2%. The result is also displayed in 
a plot of Effective Modulus vs. fiber Volume Ratio in Fig. 8.3.5. The effective modulus is 
defined as the average stress/average strain. The three-dimensional results follow closely 
to the two-dimensional solution. 
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8.4 Thick Cylinder with Circumferential Fiber Supports 

The strength of a cylinder under internal pressure can be increased by adding cir- 
cumferential fiber supports. In the present example, a three-dimensional, open-end, thick 
cylinder with four fiber supports is analyzed. The inner and outer radii of the cy Under are 
10 and 20 inches respectively, the height is 2 inches, and the radius of the fully-bonded fiber 
is 0.5 inches. By using roUer boundary conditions on the faces of symmetry, only a fifteen 
degree slice of the thick cylinder needs to be modeled. As shown in Figure 7.4.1, eighteen 
quadratic boundary elements are used to define the geometry of the model (nine-node 
boundary elements are used on both the internal and external faces of the cylinder) and 
three fiber elements are used per fiber. Three analyses are carried out with three different 
elastic moduli of fibers: 100, 500, and 1000 psi. The elastic modulus of the cylinder is 
assumed to be 100.0 psi. The Poisson ratio is 0.3 for both the matrix and the fibers. The 
cylinder is subjected to an internal pressure of 100 psi. 

Two analyses are carried out for each fiber type in which the functional variation over 
the boundary elements of the outer surface assumed a quadratic variation in one case 
and a quartic variation in another. These two analyses are compared with a multi-region, 
axisymmetric BEM analysis in which quadratic boundary elements are used to model the 
outer surface, the hole surfaces, and the fiber surfaces (Figure 8.4.2). 

The radial displacement along the top face of the cylinder is shown in Figure 8.4.3. 
Both the quadratic and quartic element approaches are in excellent agreement with the 
axisymmetric solution in the analysis in which the modulus of elasticity of the fibers and 
the matrix are the same. The quadratic and quartic element results deviate slightly from 
the axisymmetric solution for the fiber-matrix modulus ratios of 5 to 1 and 10 to 1. The 
results of the quartic elements are in better agreement with the axisymmetric results than 
are the quadratic element results. 

In Figure 8.4.4, the circumferential stress is shown for points along the top of the 
cylinder. The quartic element solution is in excellent agreement for all three fiber-matrix 
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ratios, while the quadratic element solution increases in error with increasing fiber modulus 
ratios. 

In Figure 8.4.5 the radial stress is shown along the top of the cylinder. Once again the 
results of the quartic elements are superior to the results of the quadratic elements. 

8.5 Cube with Multiple Fibers with Random Orientation 

In an attempt to analyze a material with a random fiber structure, cubes with multiple 
fibers oriented in random directions are studied. The cubes are of unit length and have 
four boundary elements per side (Fig. 8.5.1a). Randomly oriented fibers of variable length 
with radii of 0.05 are placed in five cubes in quantities of 5, 10, 15, 20 and 25 (Fig. 8.5.1b- 
f). Three cases of material properties are considered for each cube. The modulus of the 
composite matrix is 100 for all cases, however, the modulus of the fibers are 500, 10,000 
and 200,000 for the three cases studied. Poisson’s ratio is uniformly 0.3 throughout. Roller 
boundary conditions are employed on three adjacent sides and a uniform normal traction 
of 100 is applied to a fourth face. 

The normal end displacement at the center of the face on the side with the applied 
traction is plotted against the number of fibers in a cube for the three materials (Fig. 8.5.2). 
The displacement decreases with increasing number of fibers per cube and increasing Ef/E 
values as expected. 

8.6 A Beam with Fiber Reinforcement in Bending 

In the last example, the applicability of the present formulation to the study of the 
micromechanical behavior of the ceramic composite is apparent. The present formulation, 
however, is equally applicable to typical problems encountered by civil engineers. Rein- 
forced concrete can now be modeled exactly as a three-dimensional body and studied in 
detail for the first time. The present example considers a reinforced concrete beam. Here 
the concrete plays the role of the composite matrix and the reinforcement bars play the 
role of the fiber fiber. In Fig. 8.6.1, a 4 x 1 x 1 beam with four fibers is modeled using 
twenty-eight quadratic boundary elements. The effect of the ratio of fiber modulus to 
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matrix modulus ( Ef/E ) is studied for a range of values between 1 and 100. The Poisson 
ratio is 0.3 for both the beam and reinforced rods. 

The beam is completely fixed at one end and a downward shear traction of 100 is 
applied to the other end. The non-dimensional vertical displacement of the end obtained 
from the present analysis is shown in Fig. 8.6.2 as a function of Ef/E. The non-dimensional 
displacement is defined as the end displacement of the reinforced beam divided by the 
displacement of a homogeneous beam under similar boundary conditions. 

The end displacement obtained from the mechanics of material solution is also displayed 
in Fig. 8.6.2 in non-dimensional form. The curvature of the two plots are very similar but 
differ in magnitude. This difference is contributed to the fact that although the mechanics 
of material solution accounts for the stiffening due to the fibers, it does not include the 
effect of interaction between fibers. 

8.7 Laminated Fiber Composite 

A laminated composite fabricated from a fiber composite material is shown in Fig. 
8.7.1. The fiber composite is constructed with a single row of fully-bonded fibers oriented 
in the same direction. A two-ply laminate is then constructed from the fiber composite 
with the fibers of the two layers oriented at 90° angles. A boundary element model created 
for the study of this material is shown in Fig. 8.7.2. A small slice containing two fibers in 
each layer is used. The model consists of two regions. The outer surface of each region is 
modeled with sixteen quadratic boundary elements and each fiber contains two quadratic 
fiber elements. The interface between the two regions is assumed to be a perfect bond, 
however, the present version of the program also allows for sliding and spring connections. 

The composite structure is subjected to bi-axial tension. This is accomplished with 
normal tractions of 100 applied to two adjacent roller boundary conditions applied to the 
opposite ends. The elastic modulus of the composite matrix of both regions are assumed 
to be 100, and the moduli of the fibers very between 100 and 10,000. The Poisson ratio is 
0.3 for both the composite matrix and fibers at all times. 
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Figure 8.7.3 displays the displacement as a function of fiber moduli for a point on the 
interface at the corner of the plate adjacent to the sides with the applied traction. The 
material exhibits less displacement as the modulus is increased, as expected. 

8.8 Heat Conduction: Cube with Random Fibers 

The cube with randomly oriented fibers shown in Fig. 8.5.1 was analyzed for heat 
conduction analysis. The left end was subjected to a prescribed temperature of 100°F and 
for the right end a temperature of 0° was specified. All other surfaces were assumed to be 
insulated. Figure 8.8.1 shows the equivalent thermal conductivity of a cube for different 
fiber arrangements. In this analysis the ratio of thermal conductivities between the fiber 
and the matrix was assumed to be 100. 

8.9 Thermoelasticity: Effective Coefficient of Thermal Expansion 

The addition of fiber fibers in a material alters the thermal expansion of the material. 
The effective coefficient of thermal expansion of a composite material is dependent on many 
factors such as: The elastic and thermal properties of the individual constituents; the size, 
shape, orientation, and number of fibers; and the interface interaction between the fiber 
and matrix. Some investigators [Hopkins and Chamis, 1985] have formulated equations for 
the determination of the effective material properties of fiber composite materials, includ- 
ing the coefficient of thermal expansion. These equations axe, however, limited to specific 
types of fiber arrangements and interface connections. For general arrangements, experi- 
ments can be carried out for the prediction of effective material properties. Experiments, 
however, axe expensive and time consuming. The finite element method may be employed 
for this purpose, however, the convergence of the solution may be slow and the solution 
may be expensive and difficult to achieve. The present uncoupled thermoelastic BEM 
implementation is an ideal alternative for the prediction of effective material properties 
for general, tubular fiber composite materials. The analysis is relatively simple and cost 
efficient. 
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A boundary element model of a cube with nine fibers is shown in Fig. 8.9.1. The fibers 
are assumed to be perfectly bonded to the composite matrix. Twenty-four boundary 
elements axe used to model the outer surface of the cube and one fiber element is used to 
model each of the nine fibers. The cube is subjected to a uniform temperature increase by 
simply specifying the temperature on one face of the cube and zero flux on the other five 
sides. The cube is allowed to expand freely, however, rigid body translation is prevented. 
The user specified radius of the fibers is easily changed to simulate various void ratio, 
therefore, minimizing the effort required for re-analysis of the cube with different fiber to 
total volume ratios. 

In Fig. 8.9.2 the effective coefficient of thermal expansions in the lateral and transverse 
directions are shown as a function of the fiber to total volume ratio for a fiber composite 
material with the following material properties: 

^matrix _ 1? Q x 10 6 ps i £ fiber _ 2 g q x 10 6 psi 

^matrix = 0.24 l/ fiber = 0.24 

^matrix _ Q 5 x l 0 ~ 6 / p Q fiber _ 2 .1 X 10~ 6 /^ 

Also shown axe the predictions by Hopkins and Chamis (1985). The solutions of the two 
analysis axe in good agreement at low fiber to toted volume ratios and deviate slightly from 
one another as the fiber to total volume ratio is increased. A graph for a similar analysis 
is shown in Fig. 8.9.3 for a perfectly-bonded fiber composite material with the following 
material properties of aluminum and steel axe used. 


^matrix _ jq q x 1q 6 pgi 

qfiber _ go o x 10 6 psi 

^matrix _ q 3 

I/ fiber = 0.3 

a matrix _ j 2 q x IQ-G/p 

a fiber = 6.0 x 10 _6 /F 


This time the coefficient of thermal expansion of the fibers is less than that of the matrix, 
and therefore, the effective coefficient of thermal expansion decreases with increasing fiber 
to total volume ratio. Once again good agreement is seen between the solutions of the two 
analyses. 
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8.10 Heat Conduction: Linear Thermal Resistant Fiber Interface 

A cube with five fibers is used to demonstrate the linear thermal resistant interface 
relation between the matrix and the fiber. The boundary element mesh containing sixteen 
surface elements and five fiber elements is shown in Fig. 8.10.1. The radius of the fibers 
is 0.1. The conductivity of the fiber and matrix are both 100.0, however, the conductivity 
across the fiber-matrix interface is varied. The cube is subjected to a temperature of 100.0 
on a face parallel to the fibers and a temperature of 0.0 on the opposite face. The other 
four sides are insulated. 

A graph of the flux at the center of the zero temperature face versus the thermal con- 
ductivity of the interface is shown in Fig. 8.10.2. For very large values of the interface 
conductivity the interface connection approaches a perfect bond and the overall solution 
approaches the solution of a homogeneous cube with a conductivity of 100.0. As the in- 
terface conductivity decreases the flux decreases. For an interface conductivity of 0.0 the 
cube exhibits the solution of a cube with perfectly insulated holes. 

8.11 Thermoelasticity: Linear Spring-Thermal Resistant Fiber Interface 

The cooling of a cube containing five fibers is analyzed using the thermoelastic formu- 
lation. Shown in Fig. 8.10.1 is a boundary element model with sixteen surface elements 
and five fiber elements of radius 0.1. The fibers are assumed to be connected by a linear 
spring-thermal resistant interface. To demonstrate the effect of the interface connection, 
the fiber and the matrix are assumed to be constructed from the same material with an 
elastic modulus, a Poisson ratio, a thermal conductivity and a coefficient of thermal ex- 
pansion of 100.0, 0.3, 100.0, and 0.01, respectively. The thermal conductivity across the 
fiber-matrix interface is held constant at 100.0, however, the spring constants between the 
fiber and the matrix are varied. A roller boundary condition is applied on two opposite 
faces of the cube and the other four sides are free to expand, however, rigid body trans- 
lation is prevented. The cube is then subjected to a uniform decrease in temperature of 
100 . 0 . 
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A graph is shown in Fig. 8.11.1 of the resultant tractions at the center of the roller 
boundary condition faces versus the interface spring constants ( k n = k t ). For large spring 
constant values the interface connection approaches a perfect bond, and since the material 
properties of the fiber and matrix are equal, the overall solution approaches the solution of 
a homogeneous cube. As the spring constants decrease, the tractions at the roller boundary 
conditions decrease. 

8.12 Nonlinear Fiber-Matrix Interface: Beam in Bending 

A boundary element model of a 4x1x1 beam with four longitudinal fibers from example 
8.6 is shown in Fig. 8.6.1. The effect of the ratio of the fiber modulus to matrix modulus 
(Ef/E) is studied for a range of values between 1.0 and 100 . 0 . The Poisson ratio is 0.3 
throughout. A nonlinear frictional interface between the fiber and matrix is assumed. The 
coefficient of friction is 0.2 and the normal tangential spring coefficients are both 1000.0 
times the elastic modulus of the matrix. 

The beam is completely fixed at one end and a downward shear traction of 100.0 is ap- 
plied to the other end. The non-dimensional vertical displacement of the end obtained from 
the present analysis is shown in Fig. 8.12.1 as a function of Ef/E. The non-dimensional 
displacement is defined as the end displacement of the fiber composite beam divided by the 
displacement of a homogeneous beam subjected to the same boundary conditions. The 
plot of end displacement obtained in example 8.6 for a reinforced beam with perfectly 
bonded fiber- matrix interface connections is also displayed in Fig. 8.12.1. The beam with 
the perfectly-bonded fibers exhibits greater stiffness than the beam with the nonlinear 
fiber-matrix interface connections, as expected. 

8.13 Effect of Poisson Ratio in Fiber Composites 

The goal of the present five year project is to develop an accurate, yet practical bound- 
ary element program for the analysis of ceramic composite material. A conventional BEM 
approach is computationally expensive and mesh generation is tedious for an analysis of a 
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solid which may have hundreds of fibers. Therefore the challenge of the present work is to 
introduce certain approximations in the formulation which will produce an efficient BEM 
analysis and will not compromise accuracy. One approximation unique to this BEM im- 
plementation is the assumption that boundary quantities vary as a trigonometric function 
around the circumference of the fiber and hole. This was shown to be an accurate ap- 
proximation through the numerical applications of this section. Secondly, the steady-state 
BEM formulations derived in Chapters 2-4 assume that the Poisson ratio of the fiber is 
equal to the Poisson ratio of the matrix. (This limitation was removed in the formulations 
of Chapters 5, however, these formulations are more expensive.) This second assumption 
is investigated in this example using a two-dimensional BEM analysis. 

A boundary element model of a cube with a single fiber is shown in Fig. 8.13.1. Sixteen 
quadratic boundary elements are used to model the outer boundary and eight elements are 
used to model both the hole and fiber. The modulus of elasticity of the matrix is 10 6 and the 
Poisson ratio is 0.25. The elastic modulus of the fiber is assumed to range from 10 4 to 10 s 
and the Poisson ratio from 0.0 to 0.5. The cube is secured by roller boundary conditions on 
two adjacent sides, a traction of 10 5 is applied on a third side, and the normal displacement 
is measured at the center of the remaining face. The ratio of fiber diameter to side of cube 
is 0.2 for the volume of fiber to total volume ratio of 0.031. The displacement versus the 
variation in Poisson ratio (v F /v M ) is shown in Fig. 8.13.2 for the elastic modulus ratios 
(E f /E m ) of 0.01, 0.1, 1.0, 10.0, and 100.0. The displacement in these curves are normalized 
by the displacement obtained for the respective elastic modulus ratios when the Poisson 
ratio of the fiber and matrix are the same {v F fv M = 1.0). The displacement is observed to 
be minimal in all cases. The effect is greatest for the modulus ratio E F /E M of 1.0 at the 
extreme values of v F /v M . The maximum error, however, is less than 4.5%. 

This same problem is analyzed for the fiber diameter to side of cube ratio of 0.6 for 
a fiber volume to total volume ratio of 0.283 as shown in Fig. 8.13.3. The normalized 
displacement versus v F /v M for five modulus ratios is shown in Fig. 8.13.4. A very minute 
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change with variation of v F fv M is observed for modulus ratios of 0.01 and 100.0 which are 
almost indistinguishable from one another on the graph. The variation in Poisson ratio 
is seen to have a modest effect on the displacement values for modulus ratios of 0.1 and 
10.0. For an elastic modulus of 1.0, the effect of Poisson ratio is somewhat tolerable for 
the midrange values of v F /v M , however, it has a strong effect for extreme values of v F /v M . 

Next the effect of the Poisson ratio in a cube with multiple fibers interacting with one 
another is investigated. A cube with nine fibers (Fig. 8.13.5) is analyzed under similar 
conditions as the analyses of a cube with a single fiber. The volume of fiber to total volume 
ratio is again 0.283. In Fig. 8.13.6 the displacement versus v F fv M is shown for the multiple 
fiber composite for a modulus ratio E F /E M of 1.0. For comparison, the curve for a cube 
with a single fiber with the same modulus and fiber to total volume ratio is shown. The 
effect with the multiple fibers is just slightly less than the effect with a single fiber. 

In conclusion, we note that a difference in Poisson ratio between a fiber and a matrix 
should not be ignored if the difference in the modulus of elasticity is small and the volume of 
fiber to total volume is large. Nevertheless, fibers axe generally added to material in order 
to strengthen the material, and the elasticity modulus ratio ( E F /E M ) is usually greater 
than five. Furthermore, many composite constituents often have comparable Poisson ratio 
values. Therefore, the assumption that the Poisson ratio of the fiber must be equal to the 
Poisson ratio of the matrix is not very limiting. 

8.14 Composite with Fifty-one Fibers 

In Fig. 8.14.1 a boundary element mesh is shown for the outer surface of a rectangular 
block containing fifty small fibers and one large fiber (fiber elements not shown). The fibers 
are aligned parallel to one another in the arrangement shown in Fig. 8.14.2. Due to the 
relative size of the large central fiber, this fiber is explicitly modeled as a separate region 
in the conventional boundary element sense. The other fifty fibers are each modeled with 
one quadratic fiber element. The total number of nodes in the problem is five-hundred 
and ninety-eight. The material properties of the composite constituents axe 
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Matrix 

Fiber Type 1 

Fiber Type 2 

E m = 17.2 x 10 6 psi 

E F1 = 60.0 x 10 6 psi 

E™ = 28.0 x 10 6 psi 

i/ M = 0.24 

v F1 = 0.24 

u n = 0.24 

a M = 0.5 x 10 6 /F 

a F1 = 4.0 X 10“ 6 /F 

a n = 2.1 x 10~ 6 /F 

In Fig. 8.14.3 through 8.14.9 displacement contours 

axe shown for a model constrained 

with roller boundary conditions on 

three adjacent sides (bottom, left, and back faces) and 


subjected to a normal traction of 10 5 on the small (right side) face perpendicular to the 
fibers. In. Figs. 8.14.4, 8.14.6, 8.14.8 and 8.14.9 a uniform temperature decrease of 2400°F 
has also been applied in addition to the normal traction. Contours for displacement (of 
the top, traction free face) in the direction parallel to the normal traction is shown in 
Figs. 8.14.3 and 8.14.4. In Fig. 8.14.5 and Fig. 8.14.6 the contours are shown for the 
displacement (of the top face) in the direction perpendicular to the normal traction and 
the fibers. In Fig. 8.14.7 and Fig. 8.14.8 the contours are shown for the displacement (of 
the top face) parallel to the fibers. Fig. 8.14.9 displays a three-dimensional view of 8.14.8. 

In Fig. 8.14.10 a contour plot of displacement is shown for a model subjected to a 
normal traction on the top face perpendicular to the fibers. The displacement shown is of 
the top face in the direction parallel to the loading. The opposite face is fixed on rollers, 
and the lateral faces are traction free (rigid body translation is prevented). The fifty small 
fibers primarily strengthen the matrix uniformly, and the large fibers interacts with the 
strengthened matrix resulting in the deformation displayed in the figures. 

8.15 Transient Thermoelastic Analysis of a Cube with Fibers 

An uncoupled thermoelastic analysis of a cube with one and nine fibers is studied. A 
unit cube is discretized with six quadratic boundary elements on the surface of the cube 
and one fiber element per fiber. The radius of the fiber is 0.2 for the cube with a single 
fiber and 0.1 for the cube with nine fibers. The material properties, given in consistent 
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units axe 


£ fibei : = 100.0 

^matrix — Q 

j , fiber = 0.3 

^matrix _ 0 3 

Jk fiber = 1000.0 

jj. matrix _ j Q 

pcf** = 1.0 

^matrix _ j q 


The temperature of the cube is initially at zero degrees. On a face parallel to the length of 
the fiber, the temperature is held at 0° degrees and the face is supported by rollers. On the 
opposite face the temperature is suddenly raised to 100°. All other sides axe insulated and 
rigid body translation is prevented. The results of the flux and end displacement versus 
nondimensionalized time for the traction-free face is given in Fig. 8.15.1 through 8.15.4. 

The three-dimensional fiber composite BEM results axe compared with conventional, 
two-dimensional BEM results. In Figs. 8.15.1 and 8.15.2 the flux on the free face is 
shown for a cube with a single fiber and a cube with nine fibers, respectively. Similarly, 
Figs. 8.15.3 and 8.15.4 display the results for end displacement versus time. The fiber 
composite analysis is in good agreement with the conventional, two-dimensional results. 
Slight deviation is expected due to the 2-D assumption and the use of a somewhat crude 
3-D discretization. 

8.16 Transient Heat Conduction Analysis of a Turbine Blade 

A boundary element discretization of a turbine blade is shown in Fig. 8.16.1. Half 
symmetry is employed in this model which consists of a single region with ninety-two 
quadratic boundary elements. The model is 58.2 mm long, 13.9 mm wide, the radius of 
the base is 6.95 mm, and the tip is 1.98 mm (from the plane of symmetry) in thickness at 
the largest point. 

A transient heat conduction analysis is first carried out on a homogeneous blade with 
a conductivity of k = 0.0216. The blade initially rests in thermodynamic equilibrium at zero 
temperature. Then, a gas at a temperature of 1200° is assumed to flow over the blade 
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while a gas at a temperature of 500° surrounds the base of the blade. At the leading edge 
of the blade the film coefficient is h = 0.003395 and tapers off to h = 0.00064 at the trailing 
edge. At .the base of the blade the film coefficient of h = 0.00005 is assumed. Figures 8.16.3a 
through 8.16.7a contain contour plots describing the temperature flow through the blade 
at 1.0, 2.0, 4.0, 6.0, and 8.0 seconds. 

The blade is reanalyzed, this time with eight fibers (four fibers per side) r unni ng from 
the tip of the blade to half way through the base. The radius of the fibers is 0.15 mm 
and their cross-sectional locations are shown in Fig. 8.16.2. Five (9-noded) fiber elements 
are used to model each fiber. The conductivity of the fibers is 100 times the conductivity 
of the blade. In Figs. 8.16.3b through 8.16.7b the resulting temperature distributions are 
shown in contour plots. Overall, the higher conductivity of the fibers increases the heat 
flow through the blade. Hence, the heat at the tip of the blade is carried to the base faster 
than in the homogeneous blade. This results in a temperature distribution which is lower 
than the homogeneous blade near the tip and higher near the base. 

8.17 Heat Flow in a Cube with Multiple Fibers 

The steady-state heat flow through a unit cube with one, five, and nine fibers is ana- 
lyzed. The radii of the fibers is 0.1 and the arrangement of the fibers is shown in Figure 
8.17.1. The conductivity of the fibers is assumed to be 10 times the conductivity of the 
composite matrix. A uniform flux of 100 is applied to a side of the cube parallel to the 
fibers. The opposite side is maintained at a temperature of 0°. The remaining four sides 
are insulated. 

The resulting temperature profiles on the face subjected to the flux is shown in Figure 
8.17.2 for the three fiber arrangements. Note, the effect of the number of fibers on the 
temperature profile. As the number of fibers are increased, the overall conductivity of the 
cube is increased. Hence, the heat from the applied flux is carried away (to the face which 
is maintained at 0°) at a higher rate, resulting in lower temperatures. Also note, the local 
temperature minimums in the temperature profiles are associated with the close proximity 
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of the fibers near the flux boundary. 

8.18 Effective Conductivity in a Fiber Composite 

In problem 8.9, the effective coefficient of thermal expansion was calculated for a fiber 
composite with nine fibers. In the present example, the conductivity of heat flow in a 
similar specimen is examined. The boundary element model is shown in Figure 8.9.1. The 
fibers are assumed to be perfectly bonded to the composite matrix so that the resistance 
of the heat flow across the interface is zero. Twenty-four boundary elements are used to 
model the outer surface of the cube and one fiber element is used to model each of the nine 
fibers. The cube is subjected to a temperature of 0°F on a face perpendicular to the fibers 
and 100°F on the opposite face. The remaining four faces are insulated. The total heat 
flux through the cube is calculated and the effective conductivity in the direction parallel 
to the fibers can be found. When the boundary conditions are rearranged to produce heat 
flow perpendicular to the fibers, the effective conductivity in the transverse direction can 
be determined. , 

In figure 8.18.1, the effective conductivity in both the lateral and transverse directions 
axe shown as a function of the ratio of the fiber volume to the total volume of the specimen. 
The fiber composite has the following conductivities: 

£ matrix _ g 3 Btu / hf ft F fc fiber _ 33 q Btu / hr ft p 

The solutions are in good agreement with the approximate solutions by Hopkins and 
Chamis (1985) at low fiber to total volume ratios and deviates slightly from one another 
when the fiber to total volume ratio is increased. 

A similar analysis was carried out for a fiber composite with the following conductivi- 
ties: 

^ matrix _ 25 .0 Btu/hr ft F fc fiber = 5.8 Btu/hr ft F 

Once again, the results shown in Figure 8.18.2 are in good agreement with the solutions 
obtained with the formulas given by Hopkins and Chamis (1985). 
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8.19 Effective Modulus of Elasticity 

In problem 8.9, the effective coefficient of thermal expansion was calculated for a fiber 
composite, and in the above example (8.18) the effective conductivity was determined. In 
the present example, the effective modulus of elasticity will be calculated for a similar 
specimen with 9 perfectly bonded fibers. The effective modulus of elasticity can be calcu- 
lated in both the longitudinal and transverse directions by subjecting the cube to tension 
in the directions parallel and perpendicular to the fibers, respectively. 

In Figure 8.19.1, the effective modulus of elasticity in the longitudinal and transverse 
directions axe shown as a function of the ratio of the fiber volume to the total volume of 
the specimen. The fiber composite has the following material properties: 

^matrix _ 17 0 x 10 6 psi E fiber _ 2 g q * 1Q 6 psi 

^matrix _ 0 24 i, fiber = 0.24 

The solutions are again in good agreement with the approximate solutions by Hopkins 
and Chamis (1985) for the low fiber to total volume ratios, and deviates slightly from one 
another when the volume ratio is increased. 

A similar analysis is also carried out for the aluminum- steel composite matrix of ex- 
ample 8.9. 

^matnx _ 1Q q x 1Q 6 psi £ fiber _ gg g x 10 6 ps i 

^matrix _ g 3 i/ fiber = 0.3 

Once again, good agreement is seen between the solutions as displayed in Figure 8.19.2. 

8.20 Elastoplasticity: Cube with a Single Fiber 

The elastoplastic analysis of a unit cube with a single fiber of diameter 0.6 is analyzed 
using the fiber element approach developed in Chapter 6. In order to gage the accuracy 
of the formulation, a von Mises constitutive model was employed and the solution was 
compared to a two-dimensional plane stress boundary element analysis. The discretization 
of the fiber element model is shown in Figure 8.20.1. Six quadratic boundary elements 
are used to model the outer surface of the composite. Three fiber elements (with cylin- 
drical volume cells of diameter 0.98) were used to model the fiber. The two-dimensional 


discretization is shown in Figure 8.20.2. Three regions modeled with quadratic boundary 
elements and volume cells were used. The first region contains twelve boundary elements 
and is used to model the fiber. A second region (surrounding the fiber) contains twenty- 
four boundary elements and twelve volume cells and a third region contains twenty-four 
boundary elements and sixteen volume cells and is used to model the matrix. The elastic 
modulus of the fiber and the matrix are 60 x 10 6 psi and 15.0 x 10 6 psi., respectively, and the 
Poisson’s ratio is assumed to be 0.25 in both materials. The fiber is assumed to remain 
elastic and the matrix is assumed to have a yield strength of 10.0 x 10 4 psi. 

Two different analyses were carried out. In the first case the matrix is assumed to have 
a strain hardening parameter of 15.0 x 10 6 psi, and in the second case zero strain hardening 
is assumed. The cube is subjected to a normal traction on a side of the cube parallel to the 
fiber and a roller boundary condition is applied to the opposite face. In figure 8.20.3, the 
results of the end displacement vs. traction is shown for the two analyses. The 3D strain 
hardening iterative solution is stiffer than the 2D results. The difference in the solution is 
partly contributed to the plane stress assumption and to the fact that a finer discretization 
is used in the 2D analysis. Also, the cylindrical volume cell of the fiber element approach 
does not completely cover the entire volume of the composite matrix, and therefore, some 
of the nonlinear behavior is not accounted for. 

As a check of the iterative algorithm, the 3D fiber element analysis with strain hard- 
ening was conducted using the variable stiffness plasticity algorithm [Henry and Banerjee, 
1988b], The variable stiffness solution is in excellent agreement with the iterative solution 
(Figure 8.20.3). The analysis with zero strain hardening is also shown in Figure 8.20.3. 
Once again the 3D solution is stiffer than the 2D results, however, both solutions collapse 
close to the same load level. 

It is important to note that even though the fiber element analysis is three-dimensional, 
the computational time was slightly less than the solution time of the 2D approach. For 
problems with additional fibers, the fiber element approach will produce even greater 
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savings! 


8.21 Elastoplasticity: Cube with a Thin Fiber 

An analysis similar to the previous example is carried out for a unit cube with a fiber 
of diameter 0.2. The material properties from the previous example are used. The model 
shown in Figure 8.20.1 is employed with a cylindrical volume cell of diameter 0.98. 

A fiber composite analysis is carried out for four different hardening parameters: 15.0 x 
10 6 psi, 8.0 x 10 6 psi, 2.0 x 10 6 psi, and 0.0 psi. Figure 8.21.1 displays the end displacement vs. 
traction results for a cube subjected to a normal traction on a face parallel to the fibers 
with the opposite face constrained against normal displacement. The thin fiber strengthens 
the composite to a lesser degree than the larger fiber in the previous problem, resulting 
in a softer elastic and nonlinear solution and a lower collapse load for the zero hardening 
case. 

8.22 Elastoplasticity: Cube with Multiple Fibers 

An elastoplastic analysis of a cube with 9 fibers in a square packing arrangement is 
analyzed. The discretization of the fiber element model is shown in Figure 8.22.1. Six 
quadratic boundary elements are used to model the outer surface of the unit cube and 
one fiber element is used to model each fiber. The fibers have a diameter of 0.2 and a 
volume cylinder diameter of 0.32. The elastic modulus of the fiber and matrix is 60.0 x 10 6 
psi and 15.0 x 10 6 psi, respectively, and the Poisson’s ratio is assumed to be 0.25 in both 
materials. The fibers axe assumed to remain elastic and the matrix is assumed to have a 
yield strength of 10.0 x 10 4 psi. 

The analysis was carried out for a number of different hardening parameters. The 
cube is subjected to a normal traction on a face parallel to the fibers and a roller boundary 
condition is applied on the opposite face. In Figure 8.22.2, the end displacement vs. 
traction is shown. It is interesting to compare these results to the solution of the single 
fiber in problem 8.20 (Figure 8.20.3) which involved similar material properties, similar 
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loading, and had the same volume of fiber to total volume ratio. The analysis of the 
specimen with 9 fibers yields at a slightly lower stress level, however, exhibits a greater 
stiffness in both the elastic and nonlinear solution, than the analysis of the single fiber. 
The results of the zero hardening case indicate a lower collapse load in the present analysis 
than in example 8.20. 

8.23 Elastoplastic-Fracture: Cube with a Fiber 

The elastoplastic fracture model described in Section 6.5 is used to model the consti- 
tutive behavior of a unit cube with a fiber of diameter 0.6. The modulus of elasticity of 
the fiber and the matrix are 1000.0 psi and 100.0 psi, respectively, and the Poisson’s ratio 
is 0.3 for both materials. The parameters associated with the elastoplastic-fracture model 
axe given below in units of psi: 

Compression Zone: 

Ty = 1000.0 Tuet = 10 , 000.0 Ay = 1000.0 Autt = 2000.0 H = 1000.0 

Tension Zone: 

T y = 100.0 T utt = 5000.0 Ay = 100.0 A^t = 200.0 H = 100.0 

Two analyses are carried out. In the first analysis, the cube is subjected to a normal 
tensile traction on a face parallel to the fiber and roller boundary conditions are applied 
on the opposite face. In the second analysis, a normal compressive traction is applied in 
a similar manner. 2D and 3D approaches are used in both analyses. The discretization 
shown in Figure 8.20.1 is used for the 3D case and the discretization of Figure 8.20.2 is 
employed for the 2D approach. 

The results of both analyses are shown in Figure 8.23.1. In Figure 8.23.2, a close up 
view is shown for the tension case. The results of the 2D and 3D approaches agree quite 
well in both the tension and compression analyses. 
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8.24 Elastoplastic-Fracture: Fiber Composite under Axial Loading 

A 1x1x4 rectangular block has a single fiber at its center running through its length 
(Figure 8.24.1). The elastic and inelastic material properties of the composite of example 
8.23 axe assumed in the present analysis. The effect of the fiber radius on the deformation 
of the composite is carried out. Three cases were examined in which the radius of the fiber 
is 0 . 3 , 0.2 and 0.0 (homogeneous case). 

In Figure 8.24.2, the axial strain at the center of the specimen (inside the fiber) is shown 
for the three cases under both tension and compression loading. Significant stiffening is 
observed in both loading configurations as the radius is increased. In previous examples 
(8.20-8.23) it was observed that the presence of the fibers resulted in a significant loss 
of strength when the composite was loaded perpendicular to the fibers beyond its elastic 
limit. In the present analysis, a specimen loaded in tension yields at a very low level. 
For a specimen in compression the yield is higher and the hardening parameter is greater 
resulting in a very slight deviation from the elastic solution. Since the fiber remains elastic 
and the strain hardening is constant, the inelastic deformation is relatively linear once 
yielding occurs. 
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Fig. 8.2.1 Discretization of a Fiber in a Unit Cube Utilizing 
Quadratic Fiber Elements 



Fig. 8.2.2 Full Three-dimensional, Multi-region Discretization 
of a Fiber in a Unit Cube 
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Fig. £ .2.3 Comparison of Displacement Profiles Between the Full 
3-D Model and the Fiber Element Model for the End of 
a Cube in Tension 
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Fig. 8.2.4 Axial Stress Through the Cross Section of a Unit Cube 
in Tension with a Single Fiber in Parallel with the 
Loading 
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LATERAL DISPLACEMENT / LENGTH 


Fig. 8.2.6 Lateral Displacement Along a Side of a Cube Subjected 
to a Shear Force in the Cross-Plane of the Fiber 









Fig. 8.3.1 Arrangement of Multiple' Fibers in a Unit Cube 
Subjected to Lateral Tension 
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Fig. 8.4.3 Radial Displacement Through a Pressurized Thick 
Cylinder with Circumferential Fibers 
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Fig. 8.4.4 Circumferential Stress Through a Pressurized Thick 
Cylinder with Circumferential Fibers 
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Fig. 8.4.5 Radial Stress Through a Pressurized Thick Cylinder 
with Circumferential Fibers 
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Fig. 8.5.1 (a) Surface Discretization of a Unit Cube Used in the 

Study of Random Oriented Fibers , ( b— f ) Orientation of 
Variable Length Fibers Within Unit Cubes Containing 
5, 10, 15, 20, and 25 Fibers , respectively 
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Fig. 8.5.2 End Displacement of a Unit Cube With Random Oriented 
Fibers of Ef/E m = 5.0, 10.0 and 100.0 
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Fig. 8.6.1 Discretization of a Reinforced Beam Utilizing 

Quadratic Fiber Elements to Model the Four Fibers 
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END DISP. / HOMOGENEOUS END DISP. 



Fig. 8-6.2 Non-dimensionaliEed Vertical End Displacement of a 
Reinforced Beam in Bending vs. the Modulus of the 
Fiber Over Modulus of the Beam 



Fig. 8.7.1 Laminate-Fiber Composite 
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Fig. 8.8.1. Effective Conductivity of a Unit Cube with 
Random Oriented Fibers 


1.04 



—.67 I I 

-1.08 1.55 

Fig. 8.9.1. Discretization of a Cube with Nine Quadratic 
Fiber Elements 
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Fig. 8.9.2 


Effective Coefficient of Thermal Expansion 
in a Cube for Various Fiber Diameters 
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Fig. 8.9.3 Effective Coefficient of Thermal Expansion 
( in a Model of an Aluminum Cube with Nine 
Perfectly-bonded Steel Fibers) for Various 
Fiber Diameters 
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Fig. 8.10.1. 


Discretization of a Cube with Five Fibers 
Connected to the Matrix Through a Linear 
Constitutive Model 
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Fig. 8.10.2 Effect of a Thermal Resistant Fiber-matrix 
Interface on the Flux Through a Fiber 
Composite 
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Fig. 8.11.1 Effect of a Spring Interface Between the 

Matrix and Fiber in a Cube Subjected to a 
Uniform Temperature Decrease 
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Fig. 8.12.1 Vertical End Displacement of a Beam in 
Bending with a Nonlinear, Frictional 
Interface Between the Four Fibers and the 
Matrix 
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Fig. 8*13.1. Two-dimensional, Multi-region BEM Model of 
a Cube with a Single Fiber. 

(Fiber • Volume/Total Volume = 0.031) 
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Fig. 8.13.2. 


Effect of Poisson Ratio on the Lateral 
Displacement of a Single-fiber Composite in 
Tension for Values of E r /£^ of 0.01, 0.1, 
1.0, 10.0, and 100.0 
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Fig. 8.13.3. Two-dimensional, Multi-region BEM Model of 
a Cube with a Single Fiber 
( Fiber Volume/Total Volume = 0.283) 
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Fig. 8.13.4. 


Effect of Poisson Ratio on the Lateral 
Displacement of a Single-fiber Composite in 
Tension for values of E^/E^ of 0.01, 0.1, 1.0 
10.0, and 100.0 
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Fig. 8.14.1 Discretization of the Outer Surface of a 

Two-region Model Containing 50 Small Fibers 
(Small Fibers not shown) 
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Fig. 8.14.2 Cross-section of the Composite Displaying 
the Arrangement of One Large (Type 1) Fiber 
and Fifty Small (Type 2) Fibers 
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Fig. 8.14.3. Contours of Displacement in the Horizontal 
Direction Relative to the Page Due to a 
(Horizontal) Traction Loading 




Fig. 8.14.4 Contours of Displacement in the Horizontal 
Direction Relative to the Page Due to a 
(Horizontal) Traction and Temperature Loading 
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Fig. 8.14.5. Contours of Displacement in the Vertical 
Direction Relative to the Page Due to a 
(Horizontal) Traction Loading 



Fig. 8.14.6. Contours of Displacement in the Vertical 
Direction Relative to the Page Due to a 
(Horizontal) Traction and Temperature Loading 
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Fig. 8.14.7. Contours of Displacement in the Out of Plane 

Direction Due to a (Horizontal) Traction Loading 



Fig. 8.14.8. Contours of Displacement in the Out of Plane 
Direction Due to a (Horizontal) Traction and 
Temperature Loading 
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Fig. 8.14.10. Contours of Displacement in the Out-of-Plane 

Direction Due to a Traction in the Out-of-Plane 
Direction (Symmetric Boundary Conditions) 
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Fig. 8.15.3 Flux vs. Time: Transient Thermoelastic Analysis 

of a Cube with Nine Fibers 
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Fig. 8.15.4 End Displacement vs. Time: Transient Thermoelastic 

Analysis of a Cube with Nine Fibers 
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DISCRETIZATION GF A TURBINE BLADE 
(HALF SYMMETRY) 


Fig. 8 .16.1 Discretization of a Turbine Blade 
( Fiber Elements Not Shown) 
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Fig. 8.16.2 Cross-Sectional View of the Turbine Blade 

Displaying the Location of the Eight Fibers 
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transient thermal analysis of a turbine blade 


Fig. 8.16.3 Transient Heat Conduction Analysis of a 
Turbine Blade at 1 Second 
(Contour Plots of Temperature) 



Fig. 8.16.4 Transient Heat Conduction Analysis of a Turbine 
Blade at 2 Seconds 
(Contour Plots of Temperature) 
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Fig. 8.17.1 Arrangement of Multiple Fibers in a Unit Cube 
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Fig. 8.17.2 Temperature Profile in a Cube for Multiple 
Fiber Arrangements 
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Fig. 8.18.1 Effective Conductivity in a Cube for Various 
Fiber Diameters 
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Fig. 8,18.2 Effective Conductivity in a Cube for Various 
Fiber Diameters 
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Fig. 8.19.1 Effective Modulus of Elasticity in a Cube for 
Various Fiber Diameters 
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Fig. 8.19.2 Effective Modulus of Elasticity in a Cube for 
Various Fiber Diameters 
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Fig. 8.20.1 Discretization of a Fiber in a Unit Cube Using 
Fiber Elements (includes a cylindrical volume 
cell not shown) 
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Fig. 8.20.3 Load vs. Displacement Response in a Fiber Composite 
With and Without Strain Hardening 
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Fig. 8.21.1 Load vs. Displacement Response in a Fiber Composite 
for Various Strain Hardening Parameters 
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ig. 8.23.1 Elastoplastic-Fracture Load vs. Displacement Behavior 

in a Fiber Composite. Response under a Tensile Load and 
a Compression Load are both shown. 
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Fig. 8.23.2 Close-up of the Elastoplastic-Fracture Load vs. Displacement 
Response of the Tensile Load of Fig. 8.23.1 
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Fig. 8.24.1 Discretization of a Fiber in a Rectangular Block using 
Fiber Elements {includes a cylindrical volume cell not 
shown ) 


ELRSTIC SOLUTION 
TENSION 
COMPRESSION 




AXIAL STRAIN (ABSOLUTE VALUE) 


Fig. 8.24.2 Traction vs. Axial Strain :(in the Fiber). 

Elastoplastic-Fracture with Strain Hardening Behavior 
in Fiber Composites for Various Fiber Radii 




9. CURRENT ACHIEVEMENT 

9.1 Fourth Year Development 

The main focus of the work of this past year has been the development and implemen- 
tation of a nonlinear formulation that allows yielding of the composite matrix about the 
fibers. This and other important developments of the past year are outlined below. 

9.1.1 General Development 

New Quartic Boundary Elements for the Outer Boundary of the Composite Matrix : 
The fiber composite imposes a strain discontinuity on the fiber-matrix interface. This 
in turn results in a displacement variation which, although continuous, is quite irregular. 
To help capture the irregular variation in displacement and in the resulting tractions, 
quartic boundary elements were introduced to model the functional variation on the outer 
boundary (geometry remains quadratic). A quartic element has more degrees of freedom 
than a quadratic element. However, the computational time of an analysis using quartic 
elements is less than an analysis using quadratic elements when the total number of degrees 
of freedom are the same since fewer elements are required in the quartic element model 
and therefore less time is required for their integration. Furthermore, for the same number 
of degrees of freedom, an analysis using the quartic boundary elements is, in general, more 
accurate than an analysis using the quadratic elements. 

A problem employing quartic boundary elements was demonstrated in example 8.4 of 
Chapter 8. 

9.1.2 Nonlinear Material Behavior of Ceramic Composites 

Development of the Integral Equation for Displacement Rates : A formulation of an 

integral equation for displacement rates in a fiber composite structure undergoing plas- 
tic deformation was developed. The nonlinear effects present in the composite matrix are 
brought into the integral equation via a volume integral. 
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Development of a Cylindrical Volume Cell : A cylindrical volume cell was developed 
based on a 30 node shape function. The variation of initial stress rates are quadratic in 
the axial direction of the volume cell, linear in the radial direction, and trigometric in 
the circumferential direction. The development of a 5-node circular shape function was 
required for the circumferential variation. The 30 node shape function can accurately 
approximate the distribution of initial stress rates about a fiber which were observed in 
numerical tests. 

Semi-analytical Integration of the Volume Cell : A great deal of attention was devoted 
to the analytical integration of the volume kernel multiplied by the 30 node shape function. 
To lower the cost of the numerical integration, the volume cell is analytically integrated in 
the circumferential direction. This reduces the three-dimensional volume integration to a 
two-dimensional integration which is evaluated numerically. 

Development of Integral Equations for Stress Rates : An integral equation for stress 
rates was developed in which nonlinear effects were incorporated in the equation through 
a volume integral, this volume integral is strongly singular and had to be given special 
consideration. The 30 node cylindrical volume cell was once again employed to model the 
initial stress rates in the system. To reduce computational costs, semi-analytical integra- 
tion was once again carried out in the circumferential direction. 

Development of a Stress Rate Calculation on the Fiber-Matrix Interface : All integrals 
(including surface integrals) of the stress rate equation are strongly singular and are difficult 
to integrate numerically, for points on the fiber-matrix interface. For this reason, a special 
stress rate calculation was developed for the volume cell nodes which lie on the fiber-matrix 
interface. 

Implementation of an Elastoplastic-Fracture Model for Ceramic Composites : Anelasto- 
plastic-fracture model was implemented for the nonlinear constitutive modeling of ceramic 
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composites. An important feature of this model is its ability to sustain stress in compres- 
sion many times greater than in tension. The model also includes strain hardening. 

Implementation of a Nonlinear Algorithm : In order to solve the nonlinear system, a 
nonlinear algorithm was implemented based on an iterative algorithm employed in the 
computer code BEST3D. For fast convergence, the iterative algorithm has an acceleration 
scheme in which values of initial stress rates are estimated for future time steps based on 
the past history. 

9.1.3 Maintenance and Testing 

Considerable effort has been devoted to maintaining the computer program and im- 
proving the quality and reliability of the code. Verification problems from previous years 
are regularly executed to assure their continued operation. Furthermore, the present work 
was thoroughly tested and verified. 

9.2 Summary 

Significant progress has been achieved towards the goal of developing the general pur- 
pose boundary element program ‘BEST-CMS’ for the micro- and macro-mechanical studies 
of advanced ceramic composites (Table 9.1). In the first year of the contract the frame- 
work of ‘BEST-CMS’ was composed based on the advanced boundary element program 
BEST3D. At that time the formulation for elastostatic analysis of fully-bonded fibers was 
implemented in the code. All the general purpose features in BEST3D, such as the defini- 
tion of local boundary conditions and multi-region substructuring, were retained in order 
to facilitate the analysis of real problems encountered in industry. 

In the second year, the ceramic composite formulations for steady-state heat conduc- 
tion and steady-state uncoupled thermoelasticity were developed and implemented in the 
elastostatic code developed in the first year. Analytic integration of the kernel functions 
was performed about the circumference of the fiber, and therefore, the numerical integra- 
tion of the fiber (and hole) was reduced to an efficient, one-dimensional line integration 
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along its length. An efficient assembly and solution algorithm was developed for these 
formulations similar to the algorithm used in the elastostatic analysis. This fiber assem- 
bly scheme significantly reduced the number of unknowns in the system, and therefore, 
it reduced the number of equations and the work required to solve them. Also, the dis- 
placement and the heat conduction equations of the thermoelastic analysis axe integrated 
simultaneously, but solved separately for efficiency. 

Considerable effort has been focused on the development of nonlinear interface connec- 
tions between the fibers and the composite matrix. Up to this point only perfectly-bonded 
connections were considered. The ceramic composite formulations were reformulated to 
allow the implementation of a variety of sophisticated interface relations such as spring 
connections, sliding connections with smooth or frictional resistance, progressive debond- 
ing, yielding on interfaces, and temperature dependent connections. In addition to the 
fully-bonded fiber-matrix interface, constitutive relations were derived for spring connec- 
tions and slide connections with spring- coulomb frictional resistance. These models were 
implemented in a general manner to facilitate the addition of other interface models to be 
derived in the future. Since the interface relations are nonlinear and path dependent, an 
efficient nonlinear algorithm had to be derived. The nonlinear algorithm was implemented 
in both the elastostatic and the uncoupled thermoelastic analyses. 

Overall the cost of the steady-state ceramic composite analysis is just slightly more ex- 
pensive (for a moderate number of fibers) than the cost of analyzing the matrix without the 
fibers present. The additional cost is primarily attributed to the integration of the outer 
surface of the matrix for the additional nodes on the surface of the hole containing the 
fiber, and towards the expense of solving a slightly larger system. More importantly, the 
present analysis is far more efficient than the conventional three-dimensional, multi-region 
approach which requires significantly more nodes in the discretization of the fiber and hole. 
The additional nodes require additional equations which add to the expense of integration 
and solution. Furthermore, the conventional approach requires a two-dimensional numer- 
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ical integration over the surface of the fiber and hole as opposed to a one- dimensional 
integration used in the present method. The data preparation for the present method is 
less involved than the conventional multi-region approach and the location of the fibers 
can be readily changed in a re-analysis. Moreover, the code developed in the present work 
allows up to 500 (100 per GMR) fiber elements in an analysis. An ordinary multi-region 
code would require a 500 region capability in order to compete. In terms of computer 
expense and the cost of data preparation for a 500 region problem, such analyses would 
be impractical. 

In the third year, the implementation of the transient heat conduction and the transient 
thermoelastic composite analyses was carried out. These analyses are more complicated 
than their steady-state counterparts due to the presence of the convolution integrals. Since 
the transient kernels cannot be integrated analytically, a numerical integration scheme 
was developed for the integration of the fibers. Furthermore a completely new assembly 
scheme had to be developed for the analysis. This assembly scheme, will also be adopted 
for the elastodynamic formulations. Although, the transient algorithm costs more than 
its steady-state counterpart, the present implementation is very efficient relative to the 
complex nature of the problem. 

In the fourth year of development, attention was turned towards the analysis of ceramic 
composites undergoing plastic deformation. A formulation was developed and implemented 
for the yielding of the composite matrix about fibers. An elastoplastic-fracture constitutive 
law was implemented to model the nonlinear behavior of the ceramic composite. Efficient 
procedures were developed to account for the material nonlinearities and to obtain a solu- 
tion. The solution time of the iterative nonlinear analysis, varies greatly from problem to 
problem and is considerably more expensive than the elastostatic analysis. However, with- 
out the new procedures developed in this work, an elastoplastic fiber composite analysis 
would be impractical. 

Quartic boundary elements were also implemented throughout all analyses to model 
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the displacement and traction variation on the outer boundary of the composite matrix. 
The use of quartic boundary elements improves accuracy in the fiber composite analysis 
and often leads to computational savings through the reduction of boundary elements 
required in an analysis. 

The ceramic composite formulation for elastostatics, steady-state and transient heat 
conduction, and steady-state and transient uncoupled thermoelasticity has been success- 
fully implemented and has been proven to be accurate and efficient. Furthermore, the 
development and implementation of the nonlinear interface and material nonlinearities 
was successfully carried out in a general manner. Since the boundary element method has 
already been proven successful in elastodynamic and free vibration analyses, coupled with 
the success of the present work, the plan to extend the ceramic composite formulations to 
these other areas holds great potential. In fact, the boundary element method may not 
only be the best tool, but may also be the only practical tool for the analysis of large 
ceramic composite problems encountered in industry. 
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TABLE 9.1 


Current State of BEST-CMS 


Elastostatic Analysis 

Steady-state Heat Conduction Anslysis 

Steady-state Uncoupled Thermoelasticity 

Transient Heat Conduction Analysis 

Transient Uncoupled Thermoelasticity 

Nonlinear Interface between the Fibers and the Matrix 

Yielding of the Matrix about the Fibers (Elastoplastic- Fracture Model) 

The Computer Program BEST-CMS (Version 2.0) and the Associated User’s Manual 
were Last Delivered to NASA-LEWIS Research Center in the Spring of 1990 
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10. FUTURE DEVELOPMENT 


Significant progress has been achieved towards the development of the general purpose, 
ceramic composite, boundary element code ‘BEST-CMS’. Presently, the elastostatic, the 
steady-state heat conduction, and the steady-state uncoupled thermoelastic formulations 
for the analysis of fiber composite bodies with linear and frictional fiber-matrix interface 
connections have been implemented along with the transient heat conduction and the tran- 
sient thermoelastic analysis of composites with fully-bonded interfaces. During the past 
year a formulation for the analysis of ceramic composites undergoing plastic deformation 
was implemented and tested. 

In the fifth year of this contract steady-state and transient elastodynamic analyses 
and a free vibration analysis will be developed. Due to the complexity of these analyses, 
considerable effort will be required for their development. The framework for the transient 
elastodynamic analysis, however, is already in place since it closely follows the transient 
thermoelastic implementation. 

The current algorithm for nonlinear interface connections and nonlinear material anal- 
ysis have been successfully implemented and verified, however, these analyses can not be 
used simultaneously. Additional work is planned for a combined analysis in which these 
two failure modes can be studied concurrently. Also, a nonlinear thermoplastic creep model 
will be developed for an uncoupled thermoelastic nonlinear analysis. 

Due to the complex nature of the nonlinear analysis there are a number of parameters 
that must be set appropriately. Load divisions, convergence tolerances, scale factors, aspect 
ratios, etc., must be manipulated and this requires patience and an in- depth knowledge of 
the program in order to obtain a correct solution to a nonlinear problem. With additional 
development, the algorithms can be refined to remove this burden from the user. The 
development of an enhanced, user friendly, algorithm is planned in which all parameters 
would be controlled internally. The user would only need to specify the level of accuracy 


(i.e., low, medium or high) that is, desired. The end result will allow analyses of very 
sophisticated modes of failure with relative ease. 

Table 10.1 summarizes the proposed future development of BEST-CMS. The next 
version of BEST-CMS and associated manuals will be delivered at the end of this year. 

‘BEST-CMS’ will provide a very precise, yet very efficient, user-friendly, design and 
analysis tool for ceramic composites exposed to severe operating conditions. BEST-CMS 
will enable an engineer to undertake rapid numerical experiments to gain insight into 
the micro- and macro- mechanical behavior of ceramic composite materials. Armed with 
this information, the disposition of fibers can be selected to optimize performance under 
inelastic, thermal and dynamic loading. s 
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TABLE 10.1 


Planned Development of BEST-CMS 

Steady-state Elastodynamic Analysis 
Transient Elastodynamic Analysis 
Free Fibration Analysis 

Inclusion of a Temperature Dependent Creep Model in the Nonlinear Analysis 
(Thermoplastic-Creep Model) 

Combine Nonlinear Material Behavior with Nonlinear Interface Conditions 

Enhancements to Make Nonlinear Analysis User Friendly 

Delivery of the Code BEST-CMS and the Associated User’s Manual 
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APPENDIX A - STEADY-STATE KERNEL FUNCTIONS 


This appendix contains the three-dimensional, steady-state kernel functions utilized in 
Chapters 2-4. 


A.l Heat Conduction 


For the temperature kernel, 



whereas, for the flux kernel, 



A.2 Elastostatics 


For the displacement kernel, 


whereas, for the traction kernel, 

* - jsirb [- (* W -) - 
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A. 3 Steady-State Thermoelasticity 


For the displacement kernel, 




G id = 0 


■ l (wbs ) !(?)] 


Gee 


= _L(T) 

4irr \k J 


whereas, for the traction kernel, 

r, 1 1 f3yiyjVkn k \ /««»:»* + w»A „ „.. N 

Fii = 8*rM-4" l “ j “ V J 7 J ' " 


+ (v i )(‘-2-)] 


F ie = 0 






147 


APPENDIX B - TRANSIENT KERNEL FUNCTIONS 


This appendix contains the detailed presentations of the three-dimensional, transient 
kernel functions utilized in Chapter 5. For the time-dependent problems, the following 
relationships must be used to determine the proper form of the functions required in the 
boundary element discretization. That is, 

Gap(X - 0 = G a p(X - Z, n At) for n = 1 

G^(X - 0 = G aP (X - £, nA t) - G a0 (X 1)A t) for n > 1, 

with similar expressions holding for all the remaining kernels. In the specification of these 
kernels below, the arguments (X -£,t) are assumed. 

Note that the indices 

i, j, k, l vary from 1 to 3 
a,f3 vary from 1 to 4 
0 equals 4 

Additionally, 

Xi coordinates of integration point 

Zi coordinates of field point 

yi = Xi-Zi r 2 = yiyi. 

For the displacement kernel, 

Gi$ = 0 
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whereas, for the traction kernel 


* - ss*rb[ - ( SH !P !a ) - (i - » 

6 e = 0 

*"i=*[( a r )«4 

In the above, 
r) = rf(ct) 1 > 2 
c ~ k/pc e 

er f( z ) = fo e ~ x 'dx = 1 - erfc(z) 

h l(n) = erf (f) - ^e ~^/ 4 

^) = er/c(§) + ^ 

0s(r?) = er/c (§) 

/ 6 (r ? ) = er/c(§) + ^ 

/7(«?) = er/c (§) + 2^21 

f&(ri) = l-h 1 (r)). 


For the interior stress boundary kernels, 

tp _ 2/i^ £ 3G 0i , ( dGpi , 3G#, or ^ 

Epij - r^ 6ij dh v ^ 36 J 

_ 2/it/ 5F/3! ..fdFpi dFpj \ 

" 1 - 2t/ ,J 36 + P V 36 36 ) P 6iiFpe 


where 


dGij l l I" / Zymyk _ &jkyi bikyj \ 

36 167rr 2 /i(l — v) [ \ r 3 r r ) 


+ (~~) ( 3 " 4,/ ) 


3G e 

36 


1 - d? (sccfso) I W *«» - W) -s: 
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dFjj 1 1 r / I5yiyjy k yini 3 yiyjn k Z6j k ymni 

d£ k 8 tt r 3 (1 — v) [ \ r 4 r 2 r 2 

Mijyjywi\ (ZSaykyini c , 3 y^n,- 

J ~ ^ + — «* n J 

+ 5^- «,.*)<! -*)] 

dF ej 1 ( P \\{Zl IjVkVini yj n k 6 jk yini\ t fVknj\ , ,_ N 

= Ur^J [\— r ~ ~v~) f6(n) ~ ( — j m 

and the prime, ', represents a derivative with respect to r?. Thus, 

„ dfejn) 

fe dr) ' 


150 


APPENDIX C - LIST OF SYMBOLS 

A short list of notation used in this report is given below. Bold print is used to identify 
vectors and matrices. 


A b ,B b , C b 
K a ,B a ,C° 
CiAt) 
di 

Df jk i,D% kl 


Cij, F ij 


l.ep 

K ij 

L p (v UV2) 
N^V) 
M a (0 ) 


Ui 


x,X 




'» 3 




Boimdary system matrices in assembled form 
Matrices of the stress rate equations 
A tensor dependent on location of the field point £ 

Difference in displacement (and/or temperature) across the matrix-insert interface 
The elastic and the elastoplastic constitutive matrix relations, respectively 
Kernel functions 

Nonlinear interface constitutive relationship 

Two-dimensional shape function 

One-dimensional shape function 

Circular shape function 

Time 

Traction 

Displacement 

Refers to global coordinates of an integration point 
System vector of unknown boundary quantities 
System vector of known boundary quantities 
Strain 

Refers to local coordinates of an integration point 

refers to coordinates of a field point 

Stress 

Initial stress rates (corrective stress rates) 
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Subscript 


i, j, k Indicia! notation 

i,j, k = 1 for heat conduction analysis (Chapter 2) 
i,j,k = 1 , 2,3 for elastostatic analysis (Chapter 2) and for elastoplastic analysis (Chap- 
ter 6) 

i, j, k = 1,2, 3,4 for thermoelastic analysis (Chapters 3 and 4) 
a,fi = 4 for heat conduction (Chapter 5) 
a,/? = 1,2, 3, 4 for thermoelasticity (Chapter 5) 


Superscript 

H refers to quantities on the surface of the hole (left for fiber) in the ceramic matrix 
F, F2 refers to quantities on the surface of the fiber 

O refers to quantities on the outer surface of the ceramic matrix 
. (Incremental) time rate of change 
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